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ABSTRACT 

Stellar feedback plays a key role in galaxy formation by regulating star formation, driving interstellar tur- 
bulence and generating galactic scale outflows. Although modern simulations of galaxy formation can resolve 
scales of ~ 10- lOOpc, star formation and feedback operate on smaller, "subgrid" scales. Great care should 
therefore be taken in order to properly account for the effect of feedback on global galaxy evolution. We in- 
vestigate the momentum and energy budget of feedback during different stages of stellar evolution, and study 
its impact on the interstellar medium using simulations of local star forming regions and galactic disks at the 
resolution affordable in modern cosmological zoom-in simulations. In particular, we present a novel subgrid 
model for the momentum injection due to radiation pressure and stellar winds from massive stars during early, 
pre-supernova evolutionary stages of young star clusters. This model is local and straightforward to implement 
in existing hydro codes without the need for radiative transfer. Early injection of momentum acts to clear out 
dense gas in star forming regions, hence limiting star formation. The reduced gas density mitigates radiative 
losses of thermal feedback energy from subsequent supernova explosions, leading to an increased overall effi- 
ciency of stellar feedback. The detailed impact of stellar feedback depends sensitively on the implementation 
and choice of parameters. Somewhat encouragingly, we find that implementations in which feedback is effi- 
cient lead to approximate self-regulation of global star formation efficiency. We compare simulation results 
using our feedback implementation to other phenomenological feedback methods, where thermal feedback 
energy is allowed to dissipate over time scales longer than the formal gas cooling time. We find that simula- 
tions with maximal momentum injection suppress star formation to a similar degree as is found in simulations 
adopting adiabatic thermal feedback. However, different feedback schemes are found to produce significant 
differences in the density and thermodynamic structure of the interstellar medium, and are hence expected to 
have a qualitatively different impact on galaxy evolution. 

Subject headings: galaxies: feedback - galaxies: ISM - methods: numerical 



1. INTRODUCTION 

Galaxy formation remains one of the most important, un- 
solved problems in modern astrophysics. In large part, this 
is because galaxy evolution depends on small-scale star for- 
mation and feedback processes, which are still poorly under- 
stood. For instance, it is now well established that stellar 
feedback from young massive stars can significantly affect the 
ISM by regulating star f ormation ( Mac Low & Klessen 2004; 
McKee & Ostriker 2007, and references therein), driving tur- 



bulence QKim et al ||2001 
Joung & Mac Low|2006f 



de Avillez & Breitschwerdt 2004; 



Agertz et al. 2009a; Tamburro et al. 



2009) and generating galac tic sca le outflows ( |Martin| |1999, 
2 005] IQppenheimer & pave|2006J ). 

One of the most salient and long-standing problems of 
galaxy formation modelling is the overprediction of baryon 
masses and concentrations of galaxies compared to expected 
values derived using dark matter halo ab undance matching 



(Conroy & Wechsler 2009 ; Guo et al. 2010 ), satellite kinemat 
ics ( jKlypin & Prada|200 9 ; More et al.|2010j), and w eak lens- 
ing ([Mandelba um et al.|2006| ), see|Behroozi et al. (20 10]) for 
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a comprehensive discussion. It is widely thought that the low 
baryon masses of galaxies are due to galactic winds driven 
by st ellar feedback at the faint end of the stellar mass func- 
tion ( |Dekel & SiIk|[T986t |Efstathiou||2000|) and by the active 
galactic nuclei (AG N) and the bright end ( (Silk &Rees| 19981 
|Benson et al.|2003j ). 

Pion e ering numerical ga l axy formatio n studies by |Katz| 
( [1992] ), |Navarro & Whrte| ( [T993) and |Katz et al.| 
demonstrated that thermal energy from SNe inefficiently cou- 
pled to the simulated ISM due to efficient radiative cooling in 
dense star forming regions. To avoid such radiative losses, an 
ad hoc delay of gas cooling in regions of recent star formation 
is often adopte d ( Gerritsen 1997 ; Thacker & Couchman 2000 , 
2001 i lStinson et al. |2006[ |Governato et al. |2007H Agertz et al'| 
2009bT Pni 1 \ jGuedes et al. ||20 1 1 1 I Scannapieco et al. ||20 1 2] )7 
which is justified by the fact that the multiphase structure of 
star forming regions with pockets of low-density hot gas is 
not resolved. While excessive radiative losses may indeed be 
partially due to resolution effects, one may also argue that 
such losses should increase with increasing resolution, as star 
forming regions can collapse to higher densities at higher res- 
olution. Indeed, simulations of supernova-driven blast waves 
with sub-parsec resolution show that most of thermal energy 
of hot gas is radiated away during the blast wave expansion 



( [Thornton et al.|1998||Cho & Kang|2008| ). It is thus necessary 
to consider other mechanisms of stellar feedback in addition 
to the energy-driven expansion of supernova (SN) bubbles. 

Observations of the giant molecular clouds (GMCs) host- 
ing young star clusters indicate that gas is often dispersed 
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Fig. 1 . — Specific luminosity from radiation (black solid line), stellar winds 
(solid blue) and supernovae type II (solid red). The data is generated using 
STARBURST99 assuming the Geneva high mass loss stellar tracks and 
solar metallicity. 



well before the first SNe explode (t < 4Myr). This can be 
partly due to the fact that natal GMCs are not gravitatio n- 
ally b ound (M ac Low & Hessen|2004||Padoan et al.|1997t |Q| 
|et al.|2004[[Kritsuk et al.|2007[|Dobbs et al.|2011| ). However, 
there is also plenty of evidence that dense star forming re- 
gions are destroyed by their HIT regions, driven by ionization 
at low cluster masses (Walch et al. 2012 ) and radiation pres- 
sure, i.e. momentum transfer from radiation emitted by young 
massive stars t o gas and dust (Murray "et al.|2005]), and stella r 
winds at high ( Matzner 2002; Krumholz & Matz ner||2009a| ). 
In the Milky Way, the majority of star formation is concen- 
trated in a few hundred massive GMCs (M c \ ~ 1O 6 M ), con- 
taining the ma j ority of the Galax y's molecular gas reservoir 
(Murray 2 01 1| ). |Fall et al.| ( [2010a| ) presented arguments favor- 
ing radiation pressure as a gas removal mechanism to explain 
the similarity in the sha pes of the mol ecular clump and stellar 
cluster mass functions. Mu rray et al. ( |2010| ) argued that radi- 
ation pressure is the dominant force m driving the expansion 
of Galactic HII regions and showed this explicitly in the case 
of several observed star forming regions. Multi- wavelength 
data of the evolved HII region around the dense R136 star 
cluster in the 30 Doradus region of LMC is also consistent 
with these arguments and implies that radiation pressure was 
the main driver of the HII regio n dynamics at radii < 75 pc 
(Lop ez et al.|20 10, see however Pellegrini et al . |20 IT) . 

The radiation from a young stellar population indeed car- 
ries a large amount of energy and momentum, as illustrated in 
FigureH] which shows the specific luminosity and mechanical 
power from st ellar winds an d SNII from a stellar population 
assuming the Kroupa (2001 ) initial ma ss function (IMF) cal - 
culated using the STARBURST99 code ( [Leitherer et al.|1999| ). 
A SNII event typically ejects 10M o at v ~ 3000km s , and 
the ste llar winds from young massive stars have a similar ve- 
locity ( Leithe rer et al.|1 999). Although the mechanical lumi- 
nosity of winds and SNII ejecta are two orders of magnitude 
smaller than the luminosity of emitted radition, their veloc- 



ity is also two orders of magnitude smaller than the speed of 
light. This makes the actual momentum injection rate of all 
sources comparable, 



Pmd ' 



C 



' PSNll ~ /?winds ' 



ech 



(1) 



where Lb i is the bolometric luminosity of a stellar population. 
As can be seen in Figure [T] the first SNIIe occur ~ 4Myr 
after the birth of stellar population, while radiation pressure 
and stellar winds operate immediately after birth. Further- 
more, the effect of radiation pressure can be significantly en- 
hanced in dense, dusty regions as UV photons absorbed by 
dust re-radiate in the infrared, increasing the momentum in- 
jection rate in proportion to the infrared optica l depth ri R , i.e. 
Pr&d ^ T\^Ljc (see, e.g., Gayley et al. 1995 for discussion 
of momentum deposition by trapped radiation). In the envi- 
ronments of massive star clusters and central r egions of star- 
bursts , values of ttr ~ 10-100 are plausible (Mu rray et al.| 
|2010| ), making momentum imparted by radiation pressure to 
the surrounding gas the dominant feedback source at early 
times it < 4Myr) in a stellar population. 

Early dispersal of gas can facilitate the survival and break- 
out of hot gas heated by SN blastwaves and dramatically in- 
crease the overall efficiency of stellar feedback. The mo- 
mentum injection due to stellar winds and radiation pres- 
sure may also be an important feedback mechanism in its 
own right. Indeed, radiation pressure has been suggested to 
play a significant role in regulating global star formation in 



galaxies and launching galactic- scale win ds ( |Haehnert|1995 
IScoville et al.||200Tl |Murray et al.||2005]). Analytic al work 



by ( |Murray et al.||20il| see also |Na5TBT Silk 2009) demon- 
strated how massive star clusters can radiatively launch large 
scale outflows, provided star formation is vigorous enough 
(^sfr ^ 0.05 M yr _1 kpc~ 2 ); an attractive property to capture 
in simulations of galaxy formation. 

Radiation pressure feedback has just recently been con- 
sidered in numerical work studying isolated galactic disk s 



( [Hopkins et~aT]|20lTa| |2012a| [Chattopadhyay et al.|[20T2] ). 
In the suite of papers by Hopkins et al., an implementation 
of radiation pressure feedback was explored using smoothed- 
particle-hydrodynamics (SPH), relying on high mass (mspH ~ 
1O 3 M ) and force resolution (~ fewpc). It was argued 
that radiation feedback has a significant effect on galax- 
ies from dwarfs to extreme starbursts, where the contribu- 
tion w as most significan t in the high surface density sys- 
tems. |Wise et al. (2012) demonstrated, using an adaptive- 
mesh-refinement (AMR) radiative transfer technique in a fully 
cosmological context, how radiation pressure in the single 
scattering regime could affect star formation rates and metal 
distributions i n a dwarf galaxies in dark matter halos of 



2.0 x 10* M . |Brook et al.| ( |2012l ) and |Stinson et al.| ( |20i2| ) 
discussed the importance of "early feedback" in their SPH 
galaxy formation simulations. These authors assume that 10% 
of the bolometric luminosity radiated by young stars is con- 
verted into thermal energy of star forming gas over a 0.8 
Myr time period, which significantly affects simulated galaxy 
properties. However, it is not clear how such a scheme relates 
to the actual processes of early feedback, which are thought 
to be momentum- rather than energy-driven. 

Cosmological zoom-in simulations of individual galaxies 
adopting a force resolution of < 50- lOOpc, while reach- 
ing z = 0, are becoming increasingly common ([Agertz et al. 
2009b] |Governato efaLl[20T0l |Guedes et aLl|2011| ). At such 
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resolution, the largest sites of star formation can be identi- 
fied in simulations directly, although their internal structure 
would not be resolved. It is hence crucial to understand how 
well we can capture the global effect of stellar feedback from 
star forming regions taking into account all plausible sources 
and mechanisms of stellar feedback at the resolution level af- 
fordable in modern cosmological simulations. 

In this paper we discuss the available energy and momen- 
tum budget from stellar winds, SNe and radiation pressure. 
The latter is implemented using a novel empirically-based 
subgrid model. Using the ada ptive-mesh-refinement (AMR) 
code RAMSES ( Teyssier 2002 ), we study the impact of these 
feedback sources in idealized simulations of star forming 
clouds and isolated disk galaxies. The simulations are per- 
formed at spatial resolution ^50-100 pc, comparable to that 
of modern state-of-the-art cosmological simulations. We in- 
vestigate how the detailed impact of stellar feedback depends 
on the implementation and choice of parameters of feedback 
schemes. We also compare a "straight-injection" approach, 
where energy and momentum is deposited directly onto the 
grid, compares to widely used phenomenological methods 
where thermal feedback energy is allowed to dissipate over 
longer time scales than expected by radiative cooling. 

The paper is organized as follows. In § [2] we discuss the 
feedback budget from SNe and stellar winds and radiation 
pressure. § [^outlines the numerical implementation of stellar 
feedback in the AMR code RAMSES. In § [4] we present ide- 
alized cloud and galactic disk simulations, and discuss how 
the different sources of feedback affect global properties of 
star formation. We conclude by summarizing our results and 
conclusions in § [5] We detail the empirically-based subgrid 
model used to compute momentum due to radiation pressure 
in Appendix A, and implementation of the second energy vari- 
able in Appendix B. 

2. STELLAR FEEDBACK AND STAR FORMATION 

2.1. Stellar feedback 

Several processes are contributing to stellar feedback, as 
stars inject energy, momentum, mass and heavy elements over 
time via SNII, SNIa, stellar winds from massive stars, radia- 
tion pressure, and secular mass loss into surrounding interstel- 
lar gas. The feedback terms we aim to quantify in this section 
are: 

Energy: E tot = £ S Nii+£sNia + ^wind 
Momentum: p tot = psm + Pwmd+Pmi (2) 
Mass loss: m tot = m SN ii + ^SNia + ^wind + ^loss 

Metals: m Z ,tot = W Z ,SNII + ^Z,SNIa + ^Z,wind + ^Z,loss • 

We choose to calculate and include the contribution of all 
feedback processes at every simulation timestep At for every 
star particle formed by our star formation recipe (see § \2.3\ . 
Feedback is thus not done instantaneously, but continuously 
in specific time periods when a given feedback process op- 
erates, taking into account the lifetime of stars of different 
masses in a stellar population. We assume that each star par- 
ticle formed in our numerical simulations represents an en- 
semble of stars with a given initial mass function (IMF). For 
st ellar masses M £ [0.1 - 100] M@, we assume the IMF form 
Kroupa| ( |200l1| r 
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Fig. 2.— Cumulative momentum from stellar winds from STARBURST99 
(black lines) compared to the subgrid approximation in Equation^ (dashed 
red line). 
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for 0.1 <M<0.5M o 
for 0.5 <M< 100 M©, 



(3) 



where A normalizes <£(M) such that total mass of stars is equal 
to the initial mass of a star particle, m*^. Note that the choice 
of IMF can significantly affect the amount of stellar feedback, 
especially the total energy and momentum output from mas- 
sive stars. For example, the IMF of Equation [3] has more than 
twice as many massive stars exploding as type II supernovae 
(assuming SNII mass range of 8-4OM ), and a Chabrier IMF 
(Chabrier 2003]) thre e times as many, com pared to the more 
bottom heavy IMF of |Kroupa etaLl ( |1993| ). 



6 The IMF suggested by |Kroupa| j200l} extends to M = 0.01 M© with a 
slope of « 0.3 ± 0.7 below M = u.Us. For the purpose of stellar feedback, 



2.1.1. Stellar winds from massive stars 

Massive stars (M > 5M ) can radiatively drive strong stel- 
lar winds from their envelopes during the first 6 Myr of stellar 
evolution, reaching terminal velocities of 1000 -3000 km s _1 
(Lam ers & CassinelE] 1 1 999| ) . The kinetic energy of these 
winds is expected to thermalize via shocks. To account for the 
energy, momentum, mass, and metal injection by such winds, 
we use calculations done with the STARBURST99 code. We 
find that the dependence of energy and momentum injection 
on metallicity can be approximated by a simple functiorQand 
we use such functional form in our simulations. Although 
the fit is approximate, its accuracy is sufficient given the un- 
certainties in the u nderlying wind models (see discussion in 
|Leitherer et al.|1992j ). 

Specifically, we approximate the cumulative energy, mo- 
mentum and mass injection, in CGS units, for a stellar popu- 
lation of age u (in Myr), birth mass m* ^ (in M ) and stellar 

accounting for the low-mass range has a negligible effect for the feedback 
energy budget presented in this paper; the total number of available SNII 
events are reduced by only ~ 6%. 

7 We adopt the Geneva high mass loss stellar tracks, and fit for the provided 
metallicities Z = 2, 1,0.4,0.2 and 0.05 Z©, assuming the IMF in Equation[3] 
We assume the feedback behavior at higher and lower metallicities to follow 
the extrapolation of our fits. 
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Fig. 3. — Estimate of the infrared optical depth tir = ^irE c i as a func- 
tion of clum p/cluster mass. We adopt k\r w 3 cm 2 g _1 (rd/100/r) 2 dSemenovl 
|et al.|2 003 1 and th e data are o bservations of molecular clumps in the Milky 
Way (blue squares Fall et al. 2010b) and young star clusters (stars) (see Ap- 
pendix [Al for furt her details). The green stars show clusters in the LMC 
<Mackey~~& Gilmore 2003) and SMC, while cyan stars show star clusters 
within Milky Way. Trie remaining star symbols show additional extragalac- 
tic star clusters, including starbursts such as M82 ( solid red) and An tennae 
(red). The data for star clusters is from lPortegies Zwart et al. ( 2010) (M82 
data from Krumholz & Matzner 2009b). Although significant scatter exists, 
massive star cluster tend to have nigh surface densities, giving rise to infrared 
optical depths tir > 10 for M c \ > 10 5 M© . 



metallicity Z* (in units of solar metallicity Z = 0.02), as 



^wind — ^*,ini^l 



.Pwind — ^*,mi.Pl 



7 \ P3 



- ems for t* < t w 



P2 



gems 



for U < t w (4) 



m w ind = ini mi In f — + 1 ] — M for t* < t w 

\m J t w 

^Z,wind = ^wind for t* <t w , 

where e w = [1.9 x 10 48 ergs M" 1 , 0.50, 0.38], p w = [1.8 x 
10 40 g cm s" 1 Mq 1 , 0.50, 0.38] and m h2 = (2.4 x 10- 2 ,4.6 x 
10~ 4 ). The wind duration is t w = 6.5 Myr. 

In Figure [2] we show an example of how the momentum 
injection for a 1O 6 M star cluster, calculated using this ap- 
proximation, compares to the STARBURST99 calculation for 
different metallicities. The momentum injection agrees quite 
well for Z > O.1Z , although we do oversimplify the time 
evolution, especially at early time (f* < 3 Myr). A similar 
conclusion holds for the wind energy injection and mass loss. 



2.1.2. Radiation pressure 
The momentum injection rate from radiation can be written 

L(t) 



as 



Pmd = (Vl+ r l2 T IR)- 



(5) 



where tir is the infrared optical depth and L(t) is the lumi- 
nosity of the stellar population. The first term describes the 
direct radiation absorption/scattering, and should in principle 
be oc [1 -exp(-Tuv)]- However, given the very large dust and 



HI opacities in the UV present in dense star forming regions, 
rji « 1 . The second term describes momentum transferred by 
infrared photons re-radiated by dust particles, and scattered 
multiple times by dust grains before they escape, where 772 is 
added to scale the fiducial value of tir (i.e., in fiducial case 

V2 = 1). 

A simple, but crude, approach to account for radiation pres- 
sure feedback would be to assume that each star particle of 
mass m* is a single star cluster with luminosity Lit) = L\ (t)m* , 
where the specific luminosity L\(t) is shown in Figure [T] and 
that the infrared optical depth tir is a constant on the order 
of ~ 1-10. The total momentum injected into the ISM at 
every time step is then simply p m $ = pAt. However, this over- 
simplifies the impact of radiation pressure, as the effect is not 
expected to be o f uniform strength in st ar clusters of differ- 
ent masses (e.g., Krumholz & Matzner 2009b). This fact is 
illustrated in Figure [3] where we estimate ttr = /^irS c i using 
observational data for cluster/clu mp masses and radii , assum- 
ing K m w 3 cm 2 g- l (T d /l00K) 2 ( [Semenov et al |2QQ3b at solar 
dust-to-gas ratios (for dust temperatures of ^ 200 K, /^r > 
5 cm 2 g _1 ). Although the scatter is significant, this rough es- 
timate illustrates that very large values of the infrared optical 
depth are plausible in massive star clusters; e.g., the observed 
densities of the star clusters in M82 allows for tir ~ 10- 100. 
In less massive star clusters (M c \ ~ 10 2 - 10 4 ), tir is of or- 
der unity and photoionization is the dominant so urce of radia- 
tive feedback (see e.g. recent numerical work by |Walch et al.j 
|2012| , although radiation pressure may be important source of 
momentum even the single scattering (ttr = 0) regime (e.g., 
|Murray et a"Ll|20TTt |Wise et al.|[20T2| ). Note that these esti- 
mates assume a homogeneous and static distribution of dense 
gas around the young star clusters, and the effective values 



of tir around young clusters are quite uncertain (e.g., Hop 



kins e t al.|20lTb{ |Kuiper et al.|2012| [Krumholz & Thompson 



2012) 



In our fiducial simulations we use a subgrid model of radi- 
ation pressure, based on conservative empiric al estimates of 



tir. Thi s approach differs from recent work by Hopkins et al. 
( |2011b| ) where attempt is made to calculate the optical depth 
directly from the density structure of the numerical simula- 
tions. The resolution of our simulations is matched to the typ- 
ical resolution of modern state-of-the-art cosmological simu- 
lations and at such resolution the density field on the scale of 
star clusters is not resolved. 

In essence, a star particle formed via the adopted star for- 
mation prescription is assumed to consist of an ensemble 
of star clusters situated in an ensemble of natal molecular 
clumps. Via the time evolution of the bolometric luminosity 
of each star cluster, calculated using STARBURST99, we ob- 
tain the momentum injection rate exerted onto each molecular 
clump. By adopting a cluster/clump mass-size relation and 
mass function compatible with observations, we then com- 
pute the total momentum injection rate p md as the integral 
over all star cluster masses represented by the star particle 
at each simulation time step. The full description of the sub- 
grid model, and the adopted fiducial parameters, is presented 
in Appendix [A| 

2.1.3. Supernovae type II 

We calculate the time at which a star of mass M ends its 
H and He burning phases, and leaves the main sequence, us- 
ing the stellar age-mas s-metallicity fit given by equation 3 in 
Raiteri et al. (1996). By inverting this equation, we obtain 
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the stellar masses exiting the main sequence at a given age 
and metallicity. At each simulation time step, At, we calcu- 
late the stellar masses A/^ and Mt*+At that bracket the stellar 
masses exiting the main sequence over the current At. If the 
masses are in range of 8-4OM , we assume they undergo 
core-collapse and end up as SNII events. The number of SNII 
events is hence given by 



NsNII - 



$(M)dM. 



(6) 



Initially, the SNII explosion energy is in the form of ki- 
netic energy of ejecta, with a typical average value of £snii = 
10 51 ergs, which is thermalized via shocks. The total thermal 
energy injected by SNII is thus 

£snii =Nsnii£snii- (7) 

The SNII ejecta also carry momentum initial momentum, 
which should be accounted for explicitly. We assume each 
supernova event imparts momentum equivalent to an ejecta 
mass m e j = 12M ejected at v e j = 3000 km s -1 , amounting to a 
total release of 

/?SNII=NsNII^ejVej (8) 

per time step. We find that the values for the amount of energy 
and momentum injected over « 40Myr, as computed above, 
are in good agreement with the total momentum and energy 
injection computed using the STA RBURST99 code. 

Following Raiteri et a L](|1996|), we ado pt the following fits 
to the results of Woosley & Weaver (1995) for the total ejected 
mass (m e j), as well as the ejected mass in iron and oxygen (mp e 
and mo), as a function of stellar mass M (in Af©): 



m ej =0.77M L06 

m Fe = 2.8x 10" 4 M L86 

m =4.6x 10" 4 M 272 



(9) 



The total and enriched amount of ejecta released at a given 
time step becomes 



™A4j,Fe,o = J m^oQOWM. (10) 



In the RAMSES implementation, we do not track separate vari- 
ables of metal species, but simply one averaged metal density 
variable. The total mass of metals returned to the ISM, ac- 
counting for the pre-existing metallicity Z* of the stellar pop- 
ulation, is 



^z,snii = (m Fe +m )(l-Z*) + m e jZ*. 



(11) 



After each feedback step the ejecta and metal mass is returned 
to the ISM, and the star particle mass is updated accordingly. 
A more sophisticated numerical treatment of chemical enrich- 
ment must ultimately include contributions from all relevant 
species, e.g . C, N, Ne, Mg, Si, Ca and S (see e.g. Wiersma 
et al.||2009| ), which we leave for a future investigation. Note 
that oxygen dominates the ejected heavy elements by mass. 

For the IMF given in Equation [3j a stellar population of 
birth mass m*^ = 1O 4 M and Z = Z produces ~ 101.4 
SNII, ejects m e j ~ 899.4Af of material and expels m Z sNii ~ 
143.8M of metals into the ISM (of which newly produced 
iron and oxygen accounts for ~ 128.4Af ). 

In addition to the SNII feedback budget discussed above, 
which can be regarded as initial injections of energy and 



momentum into the ISM, late time evolution of supernova 
remnants can in principle inject significantly more momen- 
tum. During the first ~ 10-100 years after a SNII explo- 
sion, when SN ejecta move balli stically, the adiabatic S edov- 
Taylor (S-T) stage sets in (e.g. |Ostriker & McKee 1988]), as 
the swept up inter-stellar material greatly exceeds the ejecta. 
The shock velocity is high, leading to an approximately adi- 
abatic, energy conserving evolution. After ~ 10 4 years, the 
shock wave slows down sufficiently for the cooling time of 
post- shock gas to be of the order of or less than the age of 
the remnant, and an adiabatic assumption is no longer valid. 
|Blondin et al.| ( |1998| ) calculated the transition time at which 
the cooling time equals the age of the remnant (t coo \ = £sn) to 

be « 2.9 x 10 4 E A 5 [ 11 rf^ 11 yrs, where no is the ambient density 
and £51 the thermal energy in units of 10 51 ergs. At this time, 
the momentum of the expanding shell is approximately 

p ST =M st v S t ~ 2.6 x 10 5 E^ /I7 n~ 2/17 M Q kms~ l . (12) 

Note that p^j depends very weakly on the surrounding gas 
density and linearly on E$i and may hence be ~ 5 - 15 times 
greater than the initial ejecta momentum psmi in the density 
range ^=100-0.01 cm -3 . We regard pst as an upper limit to 
what a single SN explosion can generate, as a s ubstantial por- 
tion of the energy is lost in shocks (see § |2.2.1| ), and the clas- 
sical S-T solution assumption of a perfectly intact thin shell 
expanding into a homogeneous medium is almost certainly 
a simplification. If stellar winds and radiation pressure are 
sufficiently effective in expelling gas from young star clusters 
during the first 3-4 Myrs, hot gas may simply escape the natal 
cloud via the cleared channels. A spherical model for blast- 
wave evolution is clearly incorrect in such cases. Keeping 
this in mind, a scenario of maximally efficient S-T momentum 
generation can be modelled by replacing our fiducial c hoice 
Psmi by Pst (e.g., as is done by Shetty & Ostriker 2012| ). 

2.1.4. Supernovae type la 



Following Raiteri et al. (1996), we assume that progeni- 
tors of SNIa are carbon plus oxygen white dwarfs that accrete 
mass from their binary companions. Stellar evolution theory 
predicts that the binary masses that can given rise to white 
dwarfs exceeding the Chandrasekhar limit to be in the range 
of ~ 3 - 16 M . The number of SNIa events within a star par- 
ticle, at a given simulation time with an associated time step 
At, is then 

Af S NIa= / £(M 2 )dM 2 , (13) 



where <E>(M 2 ) is the IMF of the secondary star ( Greggio & 
|Renzini|1983l|Raiteri et al.|1996| ), 



<£>(M 2 )=A 



7 



(i) », 



(14) 



where Af B is the mass of the binary, Mm = max(2Af 2 ,3Af ) 
and Af sup = M 2 + 8Af . The normalization parameter is set 
to A' = 0.24 A (see Equation [3}. Each explosion is assumed 
to release £sNia = 10 51 ergs as thermal energy, hence injecting 
a total of £sNia = NsNia^SNia at each time step. We assume 
each SNIa to be at the Chandrasekhar limit (Af ch = 1.38Af ), 
and that this is the ejected mass upon explosion leading to 

^SNIa = A^SNIa^ch- 

We allow each SNIa event to produce O.76Af of metal 
enriched material (O.13Af of 16 and O.63M of 56 Fe) 
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(Thielemann et al. 1986 ). Note that we explicitly account for 
late time mass lass of low mass stars (M < 8M Q ) until the 
point they exit the main sequence, see §2.1.5| 

For the assumed value of A', approximately 15% of all SNe 
are of type la over the lifespan of a 1M star (10 Gyr). This 
is compatible with the notion that 10-20% of the SNe rate in 
galaxies with ongoing star formati on, such as late type spirals 
(Sbc- Sd), are due to type la events ( van den Bergh & McClure 
[T994] ). 

2. 1 .5. Stellar mass loss by low mass stars 

Although low mass stars (M < 8M ) contribute a negligi- 
ble amount to the total momentum and energy budget, they 
shed a considerable amount of mass during the a symptotic gi- 
ant branch (AGB) phase o f their evolution (e.g. |Hurley et al.| 



2000). 



Kalirai et al. ( 2008 ) provides relation between the ini- 
lar mass and final mass of the remnant in the relevant 



tial ste; 
mass range: 

M fina i = (0. 109 ±0.007)Mi nitia i + 0.394 ± 0.025 Af©. (15) 

Using the average values, the fraction of mass lost from a star 
during its lifetime is 

/ioss(M initia i) = 0.891 -0.394/M initial . (16) 

Given a star particle of birth mass m*^ and age t* we calcu- 
late, at each time step At, the expelled stellar mas^Jas 

mio SS = / M/ l0SS (M)$(M)dM. (17) 



The lost stellar mass is added to the gas mass in the corre- 
sponding cell. The gas metallicity is also updated to take 
into account metals added as part of the stellar material, 
^z,ioss = Z*^ioss- The mass loss is assumed to be quiescent, 
i.e., no momentum, other than that of the natal star particle 
with respect to the ISM, or energy is released. For the IMF in 
Equation]!] a stellar population loses ~ 25% of its mass from 
stars in the mass range 0.5- 8 M during 10 Gyr of evolution. 

2.2. Feedback budget comparison 

In §[T] we stated that radiation pressure, stellar winds and 
SNe have roughly the same momentum injection rate p. This 
is shown explicitly in the left-hand panel of Figure [4] where 
we plot the time evolution of the integrated specific momen- 
tum injected into gas, i.e. j p(t)/m*dt, due to radiation, su- 
pernovae and stellar winds for Z = 0.01 Z and Z = 1 Z calcu- 
lated using formulae described above. Note that stellar winds 
and radiation pressure inject momentum into the ISM imme- 
diately after star cluster birth, while SNIIe inject momentum 
during t ~ 4.6-38Myr. The cumulative contribution of stel- 
lar winds alone dominates over SNIIe in the first ~ 14Myr 
(6Myr) for Z = 1Z (0.01 Z ). In the low metallicity case, 
~ 5 times less momentum is injected via winds into the ISM. 
As we parametrize the energy release in a similar fashion, the 
same trends are found for the shocked wind and SNe energy 
shown in the right-hand panel of Figure]?] 

At solar metallicity, the dominant source of momentum is 
radiation pressure, reaching the equivalen t total specific SN 
momentum after only 3 Myr (see Equation |A13| ), assuming a 
stellar population of mass m* = 1O 5 M . The result weakens 



8 Agert z^eTaT] {20 1 1 \ contained a typo that omitted a factor M from their 
equation 8. The actual numerical implementation was however correct. 



by a factor of ~ 2 for Z = 0.01 Z as infrared trapping be- 
comes negligible (note that we only assume photon trapping 
for t < t c \ = 3 Myr during which cluster stars are assumed to 
be fully embedded in their natal gas clump). The non-linear 
behaviour of the strength of radiation pressure with the mass 
of the stellar population is evident, as shown by comparing 
results f or = 1O 5 M and 1O 6 M , where our model (via 
Equation A13| ) predicts ttr « 26 for the latter. This illustrates 
how radiation pressure can be an important, and even domi- 
nant, source of feedback in dense gas associated with young 
massive star clusters, as it operates at early times before the 
first SNIIe explode. 

Recall that the wind and SNe momenta in Figure [4] refer to 
the initial ejecta momentum and not any late stage momentum 
generated by an expanding bubble. The momentum expected 
from the ideal adiabatic Sedov-Taylor phase (Equation 12]) is 
greater than radiation pressure momentum even in the case of 
a supermassive (M c \ = 1O 6 M ) star cluster. However, as we 
argued above, it is not clear whether the S-T solution is ap- 
plicable in the highly inhomogeneous density field of GMCs, 
especially if gas around young star clusters is partially cleared 
by early feedback. 

2.2. 1 . Thermal energy of shocked wind and SNII ejecta 

The fate of thermal energy of gas due to thermalized wind 
and SNII ejecta can b e illu strated as follows. Following 
Sutherlan d & Dopita] ( |1993| ), we define the cooling time 
scale as t coo \ = U /A net , where the thermal energy density 
U = 3nkT/2 and the net cooling function A net = n e ntA^. Here 
n = n e +rii, where n e and ni are the number density of electrons 
and ions respectively, and An is the normalized cooling rate 
in units of ergs s -1 cm 3 . 

In a fully ionized primordial plasma at T = 10 6 K, the nor- 
malized cooling rate is log(A N ) « -23.2 and « -22 for gas 
at Z = 1 Z , see Figure [5] where we plot the cooling function 
used by RAMSES in the absence of a UV background. In the 
latter case, the cooling time is 



^cool ~ 10 



100 cm" 



years, 



(18) 



and roughly 20 times greater for a pristine plasma. Clearly, 
the cooling time is very short at average densities relevant for 
GMCs, and hot gas is quickly radiated away, unless a strong 
local heating source can maintain it. 

The criterion for heating to dominate over cooling can be 
written as 

4An < P*T, (19) 
where is the mass density of stars and T is the specific 



Ceverino & 



heating rate o f the gas in units of ergs s _i M 
Klypin (2009) argued that a very large fraction of a cell's mass 
must be converted into stars for feedback heating to overcome 
the radiative cooling for gas at T = 10 4 K, where the cool- 
ing function peaks. Even at low densities, ~ 0.1 cm -3 , the 
stellar- to-gas mass fraction must be above unity. Indeed, by 
inserting typical values for cooling and heating (see Figure[T]), 
and scaling to a typical numerical resolution of Ax = 40 pc, 
Equation [l9]relation can be written as 



< 0.535 



I 1O 4 M 



/40pcV 
{~Ax-J 



10 



a ;(io 34 )' 



(20) 



where the cooling function A is in units of ergs s 1 cm 3 and 
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Fig. 4. — Left panel: Integrated specific momentum due to radiation pressure (black line), supernovae type II (blue l ines) and stellar winds (red lines) for 
stellar population of Z = 1Z© (solid lines) and Z = 0.01 Z© (dashed lines) as a function of time. A s disc ussed in § |2.1.2| the momentum transferred to gas via 
radiation pressure depends non-linearly on star population mass. Therefore, we show p m d (Equation IXl3} for ra* = 10 3 M© and 10 6 M©. For ra* = 10 5 M©, the 
total injected p m d equals the direct SNe momentum injection for Z = 1 Z©, and momentum injected by stellar winds is ~ 30% of this value. However, radiation 
pressure and stellar winds have deposited essentially all their momentum before the first SNe explode (here t ~ 4.6 Myr). This is always the case, regardless 
of metallicity. At Z = 0.01 Z©, supernova momentum dominates, as radiation pressure is weakened due to inefficient photon trapping (p m d ~ 0.6psNii after 40 
Myr). The contribution from stellar winds also becomes negligible. For ra* = 10 6 M©, the radiation pressure model predicts tjr ~ 26 (for Z = 1Z©), making 
the final momentum contribution from radiation almost an order of magnitude greater than all other feedback sources. Right panel: the integrated specific energy 
injected by the shocked SNe and wind ejecta. The lines correspond to the same metallicities as in the left panel. 
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verted into stars. This is an order o f magnitude greater than 
what is observ ed in massive GMCs (Eva ns et al.|2009[ |Mur 



io 4 io 5 io 8 
t/m [k] 

Fig. 5. — The cooling function implemented in the RAMSES code for 
Z = 0,0.01,0.1 and 1 Z© in the absence of a UV background. 



the specific heating rate T in units of ergss -1 M _1 . In a 
cell of size Ax = 40 pc, a stellar- to-gas fraction of at least 
m*/m g2LS ~ 10 is required for heating to overcome cooling 
(at T rsj few 10 4 K), which is unachievable via star formation 
alone unless at least 90% of the original cell mass was con- 



ray et al.|20 1 1 ). As argued by Ceverino & Klypin, gas cooling 
rates drop by orders of magnitude at lower gas temperatures, 
making it possible for thermal feedback to maintain greater 
pressure gradients between dense star forming regions and the 
ambient ISM. This leads to expansion of the star forming re- 
gion that lowers the average density, eventually bringing the 
medium into a regime where heating can overcome cooling. 

The estimates made above are subject to many caveats. 
While relevant to understand the fate of thermal energy in- 
jected into gas in galaxy formation simulations, the real ISM 
is multiphase and highly inhomogeneous on the scale of the 
resolution elements of such simulations. This means that 
pockets of tenous hot gas may exist within dense gas in a sim- 
ulation cell, but it also means that estimates of the cooling 
time are optimistic as they need to include a clumping fac- 
tor that is expected to be significant in star forming regions. 
However, it is unclear how efficiently t hermal energy should 
couple to the ISM in realistic settings; |Cho & Ka ng (2008 ) 
demonstrated, using high resolution simulations of SNe ex- 
plosion in pre-existing wind-blown bubbles, that less than 
~ 10% of the shocked thermal energy could be converted into 
kinetic energy, as the rest is lost in radiative shocks within the 
bubble. 

Keeping these issues in mind, the effect of gas clearing due 
to pre- SNe momentum feedback may in many situations en- 
hance the effect of feedback, which is one of the main moti- 
vations of this work. 

2.3. The star formation recipe 
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In this work we employ a fairly standard prescription star 
formation based on the star formation rate given by 



for p > p*, 



(21) 



where p g is the gas density, p* the threshold of star forma- 
tion, and £sf the star formation, or equivalently gas deple- 
tion, time. Observations indic ate that in the local universe 
^sf ~ 2Gyr (Bigiel et al. 2011), which is a manifestation of 
the fact that observed galaxies convert their gas into stars quite 
inefficiently. 

In this work we assume that ^sf = %Aff, where % = 
y / 3ir/32Gp is the local gas free-fall time and e& is the star 
formation efficiency per free-fall time. With this assumption 
Equation j2l] enforces p* ex p 1 - 5 , which is close to the observed 
projected density relation E gas ~ X ; s F r> where n ~ 1.4 (Ken 



|nicutt|1998). As noted above, 



the efficienc y of star form ation 
1% ([Krumholz & Tan 



€ff - 



is globally observed to be low, 
2007), and we discuss our adopted values of ej? in 

For now, we would like to note that the efficiency 67 star for- 
mation per free fall is usually kept fixed in galaxy formation 
simulations. However, it is likely that this is not the case in 
observations. In fact, there is ample observational and theoret- 
ical evidence for eg to depend on scal e and environment (e.g. 
Murray et al |2011[|Padoan et al.|2012| ), which will manifest 



as stochasticity of star formation efficiency. Such stochastic- 
ity can potentially have a strong impact on feedback, because 
it implies that can be high in some regions and low in oth- 
ers. The overall star formation would thus be concentrated 
in fewer star forming sites that have high star formation effi- 
ciency, even as the global star formation efficiency averaged 
over a large patch of ISM is low. We leave an investigation 
of the effects of such stochastic efficiency on the effects of 
feedback for future work (Agertz et al. in prep.), and note 
that this caveat should be kept in mind when interpreting the 

numerical result pres ented below. 

Recent work by |Gnedin et al.| ([2009]) and |Gnedin &| 



Kravtsov| ( |2011| ) relate star formation to molecular gas, hence 



Pg —> Pu 2 in Equation [2T| which can explain why metal/dust 



poor galaxies at z ~ 3, which physically should be more prone 
to H2 destruction via UV dissociation, show deviat io ns from 
the z = K-S relat ion [Gnedin & Kravtsov| ( |2010| ). |Gnedin| 
|& Kravtsov] ( |2011| ) demonstrated that the density at which 
molecular fraction reaches 50% can be approximated as 



: 25 (Z/Z©) -1 cm" 



(22) 



which we adopt in all of our simulations as the threshold for 
star formation. In addition to the density threshold we also use 
the temperature threshold by only allowing star formation to 
occur in cells of T < 10 4 K. No other conditions or thresholds 
are used. 

To ensure that the number of star particles formed during 
the course of a simulation is tractable, we sample the Equa- 
tion 121] stochastically at every fine simulation time step At. 



For a cell eligible for star formation, the number of star par- 
ticles to be formed, N, is determined using a Pois son rando m 
process (Rasera & Teyssier 2006; Dubois & Teyssier 2008) 



P(JV) = ^exp(-Ap), 



where the mean is 



X P = 



6* Ax 3 



At. 



(23) 



(24) 



Here p* is the adopted star formation rate (Equation [21]), and 
^*,min is the chosen unit mass of star particles. In this work we 
adopt m^min = W* A^max' where 77 = 0.1, and p* is taken from 

> 10 4 M^ 



Equation 22 at solar metallicity. This yields m* ?m i n ~ iu m 
for a typical resolution of Ax max = 50pc. When the Poisson 
process produces N > I star particles in a cell at a single star 
formation events, we bin these into one stellar particle of mass 
Nm*. 

3. NUMERICAL IMPLEMENTATION OF FEEDBACK 

The efficiency of stellar feedback depends not only on its 
magnitude, but also on sp ecifics of implementation in a given 
numerical code (see, e.g., Scannapie co et al.|2012| and refer- 
ences therein). In this work we are mainly interested in gaug- 
ing the impact of stellar feedback at the state-of-the-art reso- 
lution of modern galaxy formation simulat ions without resort- 
ing to ad hoc suppression of cooling (|Gemtsen|1997[|Tn acker 
& Couchman|2000l|StTnson et al.|2006t|Governato et aIT2 007 ; 
Agertz et al. 2011) or hydrodynamical decoupling of gas ele- 
ments ( Scannapieco et al. 2006; Oppenheimer & Dave 2006). 

We choose to inject energy and momentum directly into 
computational cells as follows. Over a simulation time step 
At, we calculate the thermal energy release (E tot ), as well as 
the associated mass of ejecta (m tot ) and metals (mz,tot)- These 
quantities are deposited in the 27 cells surrounding the star 
particle, although we have also carried out most of our ex- 
periments using nearest grid point approach without signifi- 
cant differences to the final result^] We explore two different 
methods to deposit momentum: 

(1) Momentum "kicks" — Over a simulation time step At, the 
momentum p tot = p^At is directly deposited isotropically in 
the 26 cells surrounding the grid cell nearest to star particle. 

(2) Non-thermal pressure — The momentum injection rate 
p to t can be thought as a non-thermal pressure corresponding 
to momentum flux through cell surface P nt = p to t/A, where 
the area A is the surface are of a cell (A = 6 Ax 2 ), or an 
arbitrary computational region, containing a young star 
particle. This pressure is calculated at every time step and is 
added to the thermal pressure, thermal to give total pressure 
^eff = ^therm+^nt that enters in the Euler equation. We describe 
this technique in detail in Appendix [B] 

The first met hod is qualitatively simil ar to what was 
considered by |Navarro & White ( |1993| ), although these 
authors compute the momentum corresponding to a fraction 
of injected SNII energy, while we specifically compute the 
momentum injection due to various specific processes that 
generate momentum. The advantage of the first imple- 
mentation method is its simplicity, as the second method 
requires minor modifications to the Riemann solver in the 



case of the MUSCL-Hancock scheme ([Toro| |1999] ) adopted 

the first 



by the RAMSES code. On the other hand, the first method 
does not explicitly affect the cell containing the feedback 
producing star particle, which will be evacuated in the case of 
a pressure-approach. We adopt the first method as our fiducial 
choice, but we present results of both implementations in ^4] 
Strong heating and/or momentum deposition in diffuse re- 
gions can lead to extremely large temperatures and veloci- 

9 One may add further sophistication to this approach by considering su- 
pernovae explosions as discrete events, hence only applying £snii + ^SNia 
wh en an integer n umber of explosions occur during over the time-step (see 
e.g.lHopkins et al.|2012b) 
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ties. To avoid this, we disallow feedback if cell temperature is 
T > 5 x 10 8 K and limit momentum feedback to deliver max- 
imum kicks of v = 1000 km s -1 . 

The effect of momentum feedback is weakened when star 
particles occur in neighboring regions, or even com putational 
cells, as momentu m cancellations will o ccur (see e.g. [Socrates] 




Hop kins et al. d201 lb ) discussed this effect 
simulations (see their figure Al and associated 
They maximized the effect of feedback by deposit- 
ing momentum isotropically from the cloud's center of mass 
found by an FOF technique. In addition, momentum was de- 
posited in a probabilistic way that ensured that each affected 
SPH particle would receive a velocity kick at the local cloud 
escape velocity. If momentum was added gradually around 
each stellar particle, akin to our current method, Hopkins et 
al. found that feedback limited star formation less efficiently 
(by factor of ~ 5 in the measured star formation histories). 
This effect should be kept in mind as an implementation un- 
certainty. 

3.1. Increasing the impact of hot gas by delayed cooling 



As we demonstrate in §4.11 momentum feedback aids in 
clearing gas out from star forming regions, and runaway heat- 
ing can occur in some regions. However, it is still not guaran- 
teed that the evolution of hot gas is a ccurate ly captured due to 
resolution effects (see discussion in § |2.2.1| ). In addition to re- 
lying solely on early momentum feedback to clear out dense 
gas, we also consider the two following methods to capture 
the maximum effect that thermal energy from SNII may have 
on their surroundings. 

The concept of allowing for an adiabatic feedback phase 
in galaxy scale simulations has been proposed by s everal au- 
thors, (see e.g. Gerritsen| |1997t |Stinson et al.||2006|), and is 
widely utilized in the community (|Governato et al.|[2007 ' 



Ager tz et al.||2011[ |Brook et aL]|2012| ). However, trie spe 



cine implementations assume the duration of this phase to 
be much longer than the ^ 10 4 years expected from ana- 



lytical arguments. Stin son et al.| ( 2006] ) proposed a scheme 
in which SNe energy is deposited in a region of size Rsp 



10 i. 74^^-0.1^-0.20 pc? where p 04 = lO- 4 /^ 1 and/) o and^o 
are the ambient pressure and density, and cooling is disabled 

} years. However, the time scale 



for u 



- 1Q6.85 /70.32 M 0.34 p-0.70 



"51 n 



1 04 



Wx corresponds to the surv ival time of the low-density cavity 
(Mc Kee & Qstriker||1977a| ), not the adiabatic phase of SNe. 
Furthermore, as supernovae energy in the Stinson et al. im- 
plementation is deli vered at every time step, as in the method 
described in §2.1.3 most of the gas in the star forming region 
will behave adiabatically for ~ 40Myr, assuming a minimum 
SNII mass of 8 M . 

Having noted that cooling suppression models typically ex- 
aggerate the effect of SNII energy they are meant to mimic, 
we consider the effects of one such model below in a subset 
of our simulations, and compare it with results of simulations 
with no delay of cooling. When a star particle forms, we as- 
sign the time variable t coo \ to a scalar in the cell containing 
the particle. This scalar field is passively advected with the 
hydro flow. At every time step, the variable is updated as 
= t[ oo[ - At. For every cell where t coo \ > 0, cooling is 
disabled. This method approxi mates the delay of c ooling im- 
plemented in SPH codes (e.g., |Stinson et al.|20 06) within the 
Eulerian hydrodynamics context]^] We explore the effects of 

10 Note however that we assign the cooling delay time t coo \ to the gas 



delayed cooling using two values of t coo \\ 10 and 40 Myr. The 
latter is, as argued above, the duration of SNe feedback for 
stars > 8M . 

3.2. Feedback energy variable 

We also investigate a scenario in which some fraction of the 
feedback energy is evolved as a separate energy variable E^, 
which is passively advected with the hydro flow and only cou- 
ples directly to the hydrodynamic flow as an effective pres- 
sure in the Euler equations. We assume that this energy dis- 
sipates over a timescale t& s , which is assumed to be longer 
than the cooling time t coo \ predicted by the cooling rates in 
dense star forming gas, see Equation [T8] This approach can 
be viewed as accounting for the effective pressure from a mul- 
tiphase medium, where local pockets of hot gas exert work on 
the surround cold phase. A lternatively, it may be viewed as 
feedback driven turbulence ( Springel 2000), although proper 
treatment of subgrid turbulence requires not only addition of 
turbulent pressure, but also significant modifications to the 
equations solved by the code i n order to accurately model tur- 
bulent stresses and dissipation ( Iapic hmo et al.|20lTj [Sch midt 
|& Fed errath 201 1 ), which is beyond the scope of this paper 

In practice, at every time step we assume that a fraction /a, 
of the total thermal feedback energy E tot is added to the feed- 
back energy and the remaining (l-ffb) enters the ther- 
mal energy of the gas. We experiment with f fb = 0.l and 0.5, 
where the lower value is mot ivated by the r adiative SN-driven 
bubble simulations of |Cho &T?a ng (2008 ). The pressure as- 
sociated with Eft enters into the sound-speed calculation, as 
well as in the Riemann solver. During the cooling step, dis- 
sipation is modelled as E^ At = E^cxpi—At/t^) in every gas 
cell. The retention of feedback energy is here rather different 
than in delayed cooling method described above; in the lat- 
ter, the cooling delay operates over fixed time only in the gas 
present in a star particle's birth region. 

Eft will only be important in local dense star forming gas, 
where most of the thermal energy is radiated away due to high 
average density. In diffuse regions the energy budget will be 
dominated by the surviving thermal energy which dissipates 
consistently on its proper cooling time scale. We assume the 
dissipation time scale to be comparable to the decay time of 
supersonic turbulence, i.e. of order of the flow crossing time 
(e.g. |Ostriker et al.|2001| ). Massive GMCs typically have sizes 
of ~ lOpc and velocity dispersions of ~ lOkms -1 , leading to 
a crossing time of t CY ~ 1 Myr. At the scale of the disk, where 
the cold gas layer thickness is an order of magnitude thicker, 
one may argue for t CY ~ 10 Myr. We hence consider feedback 
e nergy dissipatio n time in the range t& s = 1 - 10 Myr. 

Tey ssier et al.| ( |2012] ) recently demonstrated that a feedback 
scheme employing a separate energy variable, similar to what 
is described above, is quite efficient and has a significant ef- 
fect on the star formation history, and dark matter density pro- 
file, of an isolated dwarf galaxy. 

4. SIMULATIONS 

In this section we use idealized simulations of gas clouds 
and isolated galactic disks to gauge the effect of different pre- 



present in the local cell at star particle birth, while Stinson et al. ( 2006 1 assign 
f coo i to SPH particles available within the blast wave radius Ksp at every At 
for the duration of SNII explosions. The gas particles affected by delayed 
cooling at the end of the SNII phase in the Stinson et al. approach are not 
necessarily the same particles that were present at birth. Furthermore, our 
^cooi variable is allowed to mix, leading to delayed cooling in cells previously 
not associated with the young star particle's birth cell. 
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Fig. 6. — Left panel: evolution of gas density and temperature in computation cell containing star particle for the ALL (black solid line), MOMENTUM 
(blue solid line), SNMOM (magenta solid line) ENERGY simulations (red solid line), see table[^ Dashed lines indicate if Sedov-Taylor momentum was used 
instead of simply the ejecta momentum for SNIIe. Note that simulations including momentum injection from stellar wind and radiation pressure c an cl ear out 
enough gas for runaway heating before the first SNe explosion occurs (~ 4.6 Myr). The predicted transition to a runaway heating regime (Equation |20) agrees 
well with the results of simulations. Purely thermal energy feedback (ENERGY run) has almost no effect on gas density and temperature. SNe feedback with 
momentum can clear out the gas, enabling run-away heating, but considerably slower than in the ALL run. Right panel: Gas density and temperature evolution 
in cell containing star particle in the ALL simulation for different fraction of gas mass converted into stars. The black solid lines show, from right to left, stellar 
mass fraction /* = 1 , 5, 10 and 20%. The predicted transition to regime of runaway heating occurs before the first SN explosion for /* > 10%. Note that runaway 
heating and efficient clearing of gas occur even when only 1% of gas is turned into stars. 



TABLE 1 

Feedback in a cell simulations 



Simulation 


Feedback 




ALL 


E to t & Ptot, see Equation^ 




MOMENTUM 


Only momentum: p to t 




MOMENTUM_ST 


Only momentum: p tot , where psnii 


= PST 


ENERGY 


Only energy: Etot 




SN 


Only SNII energy: £snii 




SNMOM 


SNII energy and momentum: E$mh 


PSNII 


SNMOM_ST 


SNII energy and momentum: E$mh 


PST 



scriptions for stellar feedback described in the previous sec- 
tions on the local and global efficiency of star formation (the 
^SFR-^gas relation), as well as on the structural properties of 
galactic disks. A more extensive analysis of processes such as 
outflows, and study of feedback implementations in cosmo- 
logical galaxy formation simulations will be presented in fu- 
ture work (Agertz et al. in prep). The simulations considered 
here have resolution similar to the resolution of state-of-the- 
art cosmological simulations, and the results should hence be 
directly applicable to interpretation of results in cosmological 
runs. Specifically, we restrict the spatial resolution to reach 
minimum cell sizes of Ax ~ 10- lOOpc. 

4.1. Effect of feedback at the resolution scale 

It is instructive to first study the impact of feedback in a 
typical star forming computational cell. To this end, we place 
a single star particle of mass m* in a cell of size Ax = 40 pc 
within a periodic box of fixed resolution and homogeneous 
gas density p gas = lOOmncm -3 of solar metallicity. The initial 
gas temperature matters little, as it settles to a few 10 K after 
one time step. We adopt a series of (cell) stellar mass fractions 
/* = m*/m tot G {1,5, 10,20}%, consistent with observations 



of massive GMCs (Eva ns et al.|2009[|Murray et al.|201l) . In 
the case of /* = 10%, ra* = 1.6 x 1O 4 M . We are interested in 
studying the effect of feedback on the scale of individual sim- 
ulation cells, where it will be applied in actual galaxy sim- 
ulations. At such scales, the gas self-gravity is weak due to 
the softening of forces on the scale of a couple of cells, and 
is not hence not calculated properly in the actual simulations. 
For simplicity, we choose not to include self-gravity in these 
tests. 

This setup is evolved for 30 Myr using the different feed- 
back implementations shown in table[T] In all tests we employ 
momentum deposition via a non-thermal pressure term in the 
Riemann solver (method 2 in § [3]) rather than via "kicks," as 
we want to measure the impact on the central cell containing 
the star. The infrared optical depths relevant for the ab ove 
stellar mass fraction, as approximated via Equation |A11| be- 
come tir « 0.39,0.49,0.55,0.96, i.e. at most a factor of two 
boost compared to the single scattering "L/c-regime". 

In the left panel of Figure [6] we show the gas density and 
temperature evolution for the runs with /* = 10%. The evo- 
lution strongly depends on the form of feedback employed; 
while all momentum based feedback sources can evacuate the 
cell, energy-only feedback (ENERGY and SN run) has no ef- 
fect^] illustrating the common overcooling problem. If the 



11 T his result is somewhat at odds with results of |Ceverino & K lypin 
( 2009 ), who found that purely thermal feedback from winds and SNe could 
effectively over-pressurize gas of similar characteristics, leading to gas evac- 
uation. Part of the difference stems from how thermal energy is injected; 
in our simulations, energy is deposited to the gas at every time step, heat- 
ing it to ~ 10 4 K. When the new gas state is passed to the cooling routine, 
all energy is lost over one time step, bringing the dense gas back to a few 
~ 10K. Ceverino & Klypin considered thermal feedback via a heating term 
in the cooling routine, which when balanced against cooling led to a larger 
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Fig. 7. — Cloud star formation efficiency for gas metallicities of Z = 1 Zq (left) and Z = 0.01 Zq (right) in the non self-gravitating cloud simulation. Efficient 
cooling in the more enriched cloud reduces the effect of thermal energy on star formation. The opposite is true for momentum feedback, dominated by radiation 
pressure, which by itself limits e c \ to ~ 0.1. In the low metallicity cloud, gas cooling rates are lowered by more than an order of magnitude, making thermal 
energy feedback the dominant factor limiting star formation. 



initial SNe momentum, psmi, is included (SNMOM), gas is 
pushed out of the cell and the heating criterion of Equation [20 
is satisfied after ~ 20 Myr, leading to temperatures > 10 8 K. 
When all momentum sources of feedback (MOMENTUM) are 
included, the gas is effeciently evacuated from star forming 
cell, reaching n ~ 0.1 cm -3 after only ~ 10 Myr. Simulations 
adopting the more evolved Sedov-Taylor momentum (pst) for 
each SNII result in an even faster evacuation. Not surpris- 
ingly, the strongest effect on density and temperature is found 
in the ALL run, in which runaway heating set in only after 
~ 3.7 Myr due to stellar winds and radiation pressure alone. 

In the right hand panel of Figure [6] we show the evolution 
of the ALL run for different values /*. Runaway heating is 
achieved for all values of /*, even for = 1% after ~ 10 Myr. 
For /* > 10%, this occurs before the first SNIIe explode. 

We conclude that the implemented subgrid feedback pre- 
scriptions, especially early momentum injection, greatly en- 
hance the ability of feedback to disperse and heat gas in star 
forming cells, even when only ~ 1% of gas is turned into stars. 
In more realistic setups, a patch of gas will continue to form 
stars until the star cluster has destroyed its surrounding or de- 
pleted all of the gas above the star formation threshold. We 
explore these scenarios in the next section. 

4.2. Isolated cloud 

In this idealized test, a spherical cloud of dense cold gas 
(n c \ = lOOmncm -3 , T c \ = 10K) of radius r c \ = 50 pc is placed in 
pressure equilibrium with a diffuse ambient medium (ttism = 
0.1 ran cm -3 , 7} S m = 10 4 K). Star formation is then allowed 
to proceed, as described in |2.3| with eff = 10%. As we are 
interested in the behaviour of a marginally resolved ISM, we 
adopt maximum resolution of Ax = lOpc. At this resolution, 
the cloud consists of 552 cells at exactly n c \ = 100 H cm -3 , 
having a total initial gas mass of M c \ = 1.25 x 1O 6 M . 

In the following tests, we evolve the cloud with and with- 
out self-gravity, which in a very crude way can be seen as 
limiting cases of cloud virial parameter a yiv \ no self-gravity 

equilibrium temperature. 



simply means that unresolved turbulence supports the cloud 
(c^vir ^1) and vice versa. Note that we do not attempt to 
model details of star formation in giant molecular clouds, 
which requires more advanced simulation setups. Our main 
goal is simply to gauge systematic differences between dif- 
ferent feedback implementations at the resolution level that 
should be affordable in cosmological simulations in the near 
future. 

4.2.1. No self- gravity 

In Figure 17] we show evolution of star formation efficiency 
within the cloud, defined as e c \(t) =M*(t)/M c \^ where M*(t) 
is the total stellar mass formed at time t and M c i i n i is the initial 
cloud gas mass, in the simulations without self-gravity. 

For Z = 1 Z without any feedback, the cloud forms stars 
unhindered until e c \ ~ 0.75 when the cell densities fall below 
the star formation threshold. Supernovae feedback alone can 
reduce the overall efficiency to e c \ ~ 0.2. When no momentum 
from SNe is accounted for, the efficiency is somewhat larger: 
e c i ~ 0.25, and the same conclusion holds when all thermal 
energy (and no momentum) sources of stellar feedback are 
present. 

The stellar fractions differ significantly when pre-SN mo- 
mentum feedback is included. Radiation pressure alone sets 
e c i ~ 0.125, and the conversion efficiency decreases some- 
what when momentum from wind and SNe feedback is added. 
When momentum and energy deposition from all feedback 
mechanisms is included, the final efficiency approaches e c \ ~ 
0.1, although with significantly more hot gas present in com- 
parison to pure momentum feedback. The hot gas causes vig- 
orous late time expansion of the star forming region, which is 
illustrated in the time evolution of the projected density and 
temperature in Figure [5] 

The results are different in the case of low-metallicity gaj^] 
shown in the right hand-side of Figure [7] As the gas cooling 
rates are lowered, a purely energy baseafeedback scheme can 

12 We here adopt a metallicity independent star formation threshold of 
p* = 25 cm -3 to facilitate a comparison with the Z = IZ@ case. 
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Fig. 9. — Cloud star formation efficiency when self-gravity is included. 
The line types are the same as in Figure [7] Cloud contraction boosts star 
formation efficiency at early times compared to the test without self-gravity, 
but the relative impact of feedback is the same: the star formation efficiency 
is limited to much lower values when early momentum injection is included 
(lower black blue and cyan lines). 



lower the efficiency of star formation to e c \ ~ 0.1. The ef- 
fect of radiation pressure is however not much different from 
the simulation adopting Z = 1 Z , despite tjr being 100 times 
smaller (kjr oc Z). This is because ttr plays a minor role 
in both cases, as stellar masses in the local cells are small 
(m* < 1O 4 M ). 

4.2.2. With self-gravity 

In Figure [9] we show the cloud star formation efficiency for 
the self-gravitating cloud. We do not enforce hydrostatic equi- 
librium as the (unresolved) temperature profiles would imme- 
diately be erased by cooling. As the cloud now contracts, 
the global cloud star formation efficiency becomes greater by 
more than a factor of three in all simulations. However, the 
systematic trends measured in the non self-gravitating tests 



are recovered; pre-SN feedback, and specifically momentum, 
limits star formation by roughly a factor of two more effi- 
ciently than SNe feedback. 

In this setup, a stronger impact of feedback is found when 
momentum feedback is generated via a non-thermal pressure 
in the Riemann solver (method 2 in § [3j). In the left panel of 



Figure 10 we show the 'ALL" simulation adopting free-fall 
star formation efficiencies in the range eff = 0.5 - 10%. The 
10% case is here lower by a factor of two compared to mo- 
mentum feedback via "kicks". Even though we vary the star 
formation efficiency by a factor of 20, the final global conver- 
sion sta ys within e c \ ~ 7-25%, in agreement with observed 
GMCs ( [Evans et alT2009) |Murray et aT1[20TT ), compared to 
> 70% when feedback is ignored. 

This efficient self-regulation can be understood by studying 
the star formation time scale, defined as 



Pg ^cl,ini 
'SF - — - — 



(25) 



plotted in the right hand side of Figure 10 The ability for 



the cloud to contract to higher densities makes it possible to 
achieve efficient star formation regardless of initial eff; the 
star formation time-scale regulates to tsF ~ 50 Myr after which 
the cloud is destroyed by feedback. The stellar age spread in 
patch of gas is on the order of ~ 20 Myr, where the most con- 
centrated cluster of stars formed over a narrow range of a few 
Myr. This naive model is hence qualitatively in agreement 
with observed star cluster forming regions in local galaxies 
e.g. 30 Doradus, where the stars in the massive compact 
star cluster are younger than ~ 4 Myr, while the perip heral 
stars may be as old as ~ 30 Myr ( De Marchi et al.|20lT] ), see 
also conclusions by Mu rray etaL] ( 201 lj ) regarding Milky Way 
GMCs. 

It is plausible that effective self-regulation only occurs 
when simulated star forming gas cloud are resolved suf- 
ficiently for self-gravity to allow for some degree of col- 
lapse/contraction. At a cosmological resolution Ax ~ lOOpc, 
such collapse may not occur to the same degree as observed 
in the experiments here, especially as the gas is pressurized 
artificially at the scale of resolu tion to prevent spurious frag- 
mentation ( [Truelove et al.|1997] ). 
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Fig. 10. — Left panel: global cloud conversion efficiency e c \ for cloud runs with different assumed values of eff. Radiation pressure in these tests is modelled 
as a non-thermal pressure term (see text for discussion) and its inclusion lowers the over all efficiency by roughly a factor of two for eff = 10% compared to when 
momentum is deposited directly to surrounding cells. Note that although eff is varied by a factor of 20, the final stellar mass fraction only varies by a factor of two. 
The dashed line shows e c \ for a simulations with eff = 10% without feedback. Right panel: instantaneous gas consumption time scale, t$F = m*/m*. Regardless 
of the initial t$F, which simply reflects the initial conditions and choice of eff, all clouds reach t$F = 50- 100 Myr after roughly one free-fall time at which point 
stellar feedback starts to disrupt the cloud. 



4.3. Disk galaxy 

Following Hernquist (1993) and Springel (2000) (see also 
Springel et al. 2005) we create a particle distribution repre- 
senting a late type, star fo rming spiral galaxy embed ded in 
an NFW dark matter halo ( [Navarro et aTJ[T996l [T997] ). The 
halo has a concentration parameter c = 10 and virial circular 
velocity, measured at overdensity 200p cr i t , ^200 = 150kms _1 , 
which translates to a halo virial mass M200 = 1.1 x 1O 12 M . 
The total baryonic disk mass is M disk = 4.5 x 10 10 M o with 
20% in gas. The bulge-to-disk mass ratio is B/D = 0.1. We 
assume exponential profiles for the stellar and gaseous com- 
ponents and adopt a disk scale length r d = 3.6 kpc and scale 
height h = 0.1 r d for both. The bulge mass profile is that of 



4.3.1. Star formation histories 



Hernquist] ( |1990| ) with scale-length a = 0.1r d . 

We initialize the gaseous disk analytically on the AMR grid 
assuming an exponential profile. The galaxy is embedded in 
a hot (T = 10 6 K), tenuous in = 10~ 5 cm~ 3 ) gas halo enriched 
to Z = 1O~ 2 Z , while the disk has solar abundance. We con- 
duct all simulations at a maximum AMR cell resolution of 
Ax = 70 pc, typical of current state-of-the-art galaxy forma- 
tion simulations carried out to z = 0. 

We systematically vary the different sources of stellar feed- 
back operating in the simulations, and conduct additional tests 
which include thermal feedb ack via phenomenological ap- 
proaches described in § 3.1 Table [2] presents details of all 
simulations considered in the following analysis. All runs 
adop t the standard star formation prescription outlined in 
§ 2.3 and we generally adopt a star formation efficiency per 
free fall time of eff = 10%. We note that this value of effi- 
ciency is an order of magnitude larger than the average values 
derived glo bally for kiloparsec patches of gas or in ind ividual 
clouds (e.g. |Krumholz & Tan|200"7] |Bigiel et al.|20Q8j ). How- 
ever, as we show below, runs with feedback and large free-fall 
efficiencies produce normalizations of the Ken nicutt- S chmidt 
relatio n quite close to observations (see also Hopkins et al. 
|2011b| ). 



Figure 1 1 shows the star formation histories in disk simula- 
tions presented in table [2] The top left panel presents the im- 
pact of direct feedback injection, i.e. without any phenomeno- 
logical approach to thermal energy. We find the same trend in 
star formation rate as in the isolated cloud test. Simulations 
that include only thermal energy or SNe have a minor effect 
on the star formation history compared to no feedback, while 
the inclusion of momentum lowers the SFRs by up to a factor 
of three. This process is mainly due to early, pre-SN feed- 
back, especially radiation pressure. After a few orbital times 
all simulations regulate to roughly the same SFRs, although 
at different gas fractions. 

I t is instructive to compa re our results with the recent work 
by |Hopkins et al. ( 2011a ). These authors reported average 
infrared optical depths of (ttr) ~ 10-30 in their simulated 
Milky Way-like galax}[^| Our model, on the other hand, pre- 
dicts more modest average values in the range (ttr) ~ 2-6. 
The actual values of ttr in dense gas surrounding young, em- 
bedded star clusters are highly uncertain both because we do 
not know covering fraction o f absorbing dusty gas (see, e.g., 
Krumholz & Thompson |2012| ) and because dust temperatures 
used in calculations of tir are assumed to be high, T d > 100 K, 
while the optical depth can be much l ower if dust tempera - 



tures are much lower because tir oc T d 2 ([Semenov et al. 



2003). 



To understand how significantly larger values ol tir arrect 
our results, we perform two "All" simulation using fixed op- 
tical depths tir = 10 and 30. As shown in the right panel of 
Figure [TT] increasing r IR further suppresses SFR by ~ 30% 
for tir = 10, and by a factor of 2-3 for tir = 30. The latter 
case renders SFRs ~ 5 — 10 times lower than in the case of no 

13 These models marginally resolve the collapse of individual GMCs, al- 
though not their internal structure. The optical depth tir in their simulations 
steadily increases in the star forming clouds until feedback halts the grav- 
itational collapse. The reported optical depths refer to the average values, 
used in the feedback scheme, at the moment when particles are stochastically 
chosen to receive a feedback velocity "kick." 
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TABLE 2 

Galactic disk simulations of z = spiral galaxy analogue. See Equation[2]for notation. 



Run 


Description 




Direct injection runs 


nofb(eOOl) 


No feedback, e ff = 10% (1%) 




momentum 


Only momentum: p to t, eff = 10% 




energy 


Only energy: E to t, eff =10% 




prad 


Only radiation pressure: p md , eff = 10% 




SNnomom 


Only SNe energy: £snii> eff = 10% 




SN 


Only SNe energy and momentum: £snii & Psniu £ff = 10% 




all(eOOl) 


All feedback processes: .Etot & Ptot eff = 10% (1%) 




all_taulO 


All feedback processes: E tot & p tot , fixed tir = 10, eff = 10% 




all_tau30 


All feedback processes: E to t & Ptot, fixed tir = 30, eff = 10% 




Simulations adopting delayed cooling 


SN_dclO 


Only SNe energy: £sniu delayed cooling t coo \ = 40Myr, eff = 10% 




SN_dc40 


Only SNe energy: £sniu delayed cooling t coo \ = 40Myr, eff = 10% 




energy_dclO 


Only energy: E tot , delayed cooling t coo \ = lOMyr, eg- = 10% 




energy_dc40 


Only energy: E to t, delayed cooling t coo \ = 40Myr, eff = 10% 




all_dclO 


All feedback processes: E to t & Ptot, delayed cooling t coo \ = lOMyr, eff = 


10% 


all_dc40(e001) 


All feedback processes: E to t & Ptot, delayed cooling t coo \ = 40Myr, eff = 


10%(1%) 


Runs adopting a feedback energy variable 


energy_f05_tl 


Only energy: E tot , feedback energy fraction fa = 0.5, dissipation time 


s = lMyr, €ff =10% 


energy_f05_tl0 


Only energy: E to u ffb = 0.5, f dis = lOMyr, e ff = 10% 




all_f05_tl 


All feedback: E tot & p to t, ffb = 0.5, f d is = 1 Myr, e ff = 10% 




all_f05_tl0(e001) 


All feedback: E to t & Ptot, ffb = 0.5, f d is = lOMyr, e ff = 10%(1%) 




all_f01_tl 


All feedback: E tot & ptot, ffb = 0.1, *dis = 1 Myr, e ff = 10% 




all_f01_tl0 


All feedback: E tot & Ptot, ffb = 0.1, ?dis = lOMyr, e ff = 10% 




all_f01_t40 


All feedback: £ to t & Ptot, ffb = 0.1, r d is = 40 Myr, e ff = 10% 





feedback. 

In the bottom left panel we present the impact of dis- 
abling cooling in the gas surrounding newly born star parti- 
cles. The SFR in runs with t coo \ = 10 and 40 Myr is suppressed 
by amount similar to the runs with high tir values discussed 
above. A significant suppression (by a factor of two) can be 
achieved via SNe alone, provided gas cooling is disabled for 
extended periods of time, t coo \ = 40 Myr. 

The effect of treating a fraction /fb of the feedback en- 
ergy as an auxiliary energy variable that dissipates on a 
timescale t& s , longer than expected from cooling in the dense 
gas, is shown in the bottom right panel of Figure [TT] Even 
for a modest /fb = 10% dissipating over t& s = 1 Myr, SFRs can 
be affected by ~ 30%. As the energy fraction is increased 
to /ft, = 0.5, and/or dissipation occurs over longer time scale 
^dis ^ 10 Myr, we find a significant impact on the SFH s, and 
SFRs approach a steady ~ 1M yr -1 . As discussed in § 2.2.1 
up to 90% of SNe energy may be lost in ra diative shocks 



within wind-blown bubbles ( Cho & Kang 2008 ) in a few Myr. 
However, as the above simulations indicate, even this amount 
of preserved energy has a non-negligeble effect on star for- 
mation rate. We view this as an indication that some form of 
sub-grid treatment of feedback energy may be required, even 
in the presence of pre-SN feedback sources, due to the unre- 
solved ISM phases and gas motions. 

We note that these results should only be viewed as indica- 
tive, as the effect of feedback can in general depend on metal- 
licity, ISM pressure, depth of potential well, accretion rates 
etc., which we plan to explore in future work. 

4.3.2. The EsFR-5]g as relation 



per free-fall time in the presence, and absence, of feedback. 
All data points refer to quantities averaged over azimuthal 
bins of width Ar = 720pc, and are calculated from simula- 
tion snapshots in the ti me range 240 -300 Myr . Shown is also 
the THINGS data from |Bigiel eTakl (|2008 F] and t he galaxy 



scale average relation from Kennicutt ( 1998]TThe |Bigiel et al. 
(2008 ) relation is derived for kilo parsec sized patches, and is 
hence a more comparison to our simulated data. 

Without feed back, simu l ations adopting eff = 1 % are con- 
sistent with the Kennicutt (1998 ) relation. However, at high 
Sgas the adopted non-linear star formation relatio n (p* oc p 15 ) 



Figure 12 shows how the Kennicutt- Schmidt (KS) rela- 
tion is affected by the change of star formation efficiency 



over-shoots the observed, less steep relation of Bigiel et al. 
(2008J. In runs with no feedback, the normalization of the 
^sfr - ^gas relation scales linearly with the assumed value of 
eff, while in runs with feedback (the "All" model) the ampli- 
tude of the relation changes by a factor of at most two for val- 
ues of eff that differ by a factor of ten. However, data points at 
the largest values of S gas , corresponding to the galactic cen- 
ter in the analyzed simulation snapshots, are less affected by 
feedback and the difference in amplitude for runs with dif- 
ferent eff persists in these regions. We note that in runs with 
eff = 1%, the KS relation with and without feedback is similar. 

The dependency of feedback model parameters on the KS 
relation is shown in Figure 13 in which different panels show 
the effect of increasing the strength of radiation pressure, de- 
laying cooling for longer times, and increasing the contribu- 
tion/duration of feedback energy using a second energy vari- 
able. Overall, the sensitivity to the parameters is fairly weak: 
the KS relation is similar for models in which dissipation of 
SNII energy is slowed down by delay of cooling or via using 
second energy variable for £ coo i < 10 Myr or t^ s < 1 Myr, and 

14 Surface densities are corrected by a factor of 1.36 to account for helium. 
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Fig. 11. — Star formation histories for the isolated galactic disk simulation. Top left: the impact of various feedback sources, in the "straight injection" 
implementation of feedback. The strongest suppression of star formation is in simulations that include early momentum injection, especially momentum due to 
radiation pressure. Thermal energy feedback has a sub-dominant effect due to short cooling times in dense gas. Top right: the impact of increasing the infrared 
optical depth to a fixed value of tir = 10 or 30. tir = 10 reduces the SFR by ~ 30%, while boosting radiation pressure using tir = 30 suppresses the SFR by 
another factor of 2-3 compared to the fiducial "All" run, and a factor of ~ 5 - 10 compared to the case of no feedback. Bottom left: the impact of delaying cooling 
in the local gas around newborn star particles for t coo \ = 10 and 40 Myr. Note that delaying cooling for such values of t coo \ results in a similar suppression of SFR 
compared to the radiation pressure momentum injection with large values of ttr. For example, SFR w 2Mq yr _1 for t coo \ = 40 Myr, which is similar to the SFR 
for run with tir = 30 shown in the top right panel. Bottom right: the impact of assigning some fraction of the feedback energy to an energy variable Eft, that 
dissipates on longer timescales t^ s than expected from cooling in the dense gas. Even if only 10% of the energy is assumed to dissipate over ^is = 1 Myr, SFR is 
suppressed by ~ 30%. If the energy fraction is increased to /a, = 0.5, and/or dissipation occurs over longer timescale fdis > lOMyr, we find a significant impact 
on star formation, as SFR approaches a steady rate of ~ 1 Mq yr _1 . 



for models with early momentum injection with optical depth 
up to tir = 10. As parameters are dialed up to even larger val- 
ues (tir = 30, t coo \ = 40 Myr, or t^ s > 10 Myr), normalization 
of the KS relation is significantly suppressed. 

These results show that our fiducial feedback model ("All") 
at the adopted resolution level, results in star formation rates 
comparable to the runs in which cooling is delayed or SNe 
energy is dissipated on a controlled time scale. The results 
also show that normalization of the KS relation can be used 
to constrain the plausible range of values of parameters, or at 
least exclude the most extreme values. 

4.3.3. Visual comparison 



In Figure 14 we show face-on and edge-on maps at t = 
200 Myr of tEe gas surface density, mass-weighted average 
temperature within z = ±150pc of the disk, and stellar surface 
density of five of the simulations from table [2j "nofb", "all", 



"all_tau30", "all_dc40" and "all_f05_tl0". The two former 
runs are our fiducial runs with and without feedback, and the 
latter three represent efficient feedback implementations. 

In runs without feedback, dense star forming clumps of gas 
form out of spiral arms, and remain intact throughout the sim- 
ulation until star formation depletes most of their gas, or the 
clumps sink to the disk center. This run thus produces very 
massive star clusters clearly visible in the stellar surface den- 
sity map. In the "all" simulation, gas clumps do not form 
or are effectively dispersed and gas distribution in this run is 
considerably less clumpy. Consequently, massive star clusters 
are not produced, and this effect is even more pronounced in 
the three example of efficient feedback. 

All simulations feature a highly multiphase medium. Large 
holes filled with hot coronal gas at T ~ 10 6 K forms between 
the cold gas associated with the spiral arms in all simulations. 
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Fig. 12. — The impact of feedback implementation and assumed star formation efficiency on the XsFR-^gas relation. The left panel shows runs with no 
feedback, while the right panel shows runs with "All" feedback implementation (see table [2} . In both panels the two sets of points show runs with two different 
assumed star formation efficiencies eff = 1% to 10%. The points correspond to the average aisk values in azimuthal bins of wi dth Ar = 720 pc, an d are calculated 
from simulation snapshots in the time range 240-300Myr. The black solid line shows the galactic scale averaged data from Kennicutt (1998 1 and the contour 
lines the distribution of sub-kpc sized patches in the sample of nearby galaxies by ( Bigiel et al. 2008). In runs with no feedback, the normalization of the relation 
scales linearly with the assumed value of eff, while in runs with feedback the amplitude ol the relation changes by a factor of at most two for values of eg that 
differ by a factor of ten. Star formation in the central parts of the galaxy, here points with the largest values of £ gas , is not affected by feedback to the same extent 
as the rest of the disk and the difference in amplitude for runs with different eff persists in these regions. 
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Fig. 13. — The impact of different implementations of feedback on the Esfr - ^gas relation. The left panel shows the effect of varying the strength of radiation 
pressure momentum injection, the middle panel shows effect of delaying cooling around newly born star particles, and right panel shows effect of treating a 
fraction of thermal feedback energy as a separate energy variable. Data points show azimuthally averaged values adopting bin sizes of Ar = 720 pc, and are 
calculated from simulation snapshots in the time range 240-300Myr. The observational data points are described in the caption of Figure [12] Larger values 
of tir, cooling delay time t coo \, or energy dissipation time £di s , can lead to a similar suppression of normalization of the £sFR _ ^gas relation~The investigated 
feedback methods show a factor of ~ 20 spread in the normalizations of the Ssfr - ^gas relation, which shows that this relation can be a useful tool in constraining 
parameters of feedback models. 



This effect is less prominent in the simulations incorporating 
feedback, as cold gas is pushed out of star forming regions, 
resulting in a larger filling factor of cold material. This effect 
is especially apparent in the face-on temperature map of the 
"all_f05_tl0". The edge-on maps of density and temperature 
in all feedback runs show that fountains and outflows of both 
cold and warm gas (T ~ 10 4 - 10 5 K) and hot gas (T > 10 7 K) 
are present close to the disk plane. The efficient feedback runs 
all feature a more porous ISM, with prominent pockets of hot 
gas forming within spiral arms, as seen in the face-on density 



and temperature maps. 

This illustrates that specific details of feedback implemen- 
tations do matter in determining qualitative structural proper- 
ties of the ISM and even stellar distribution. We quantify the 
differences in density and temperature structure of the ISM in 
these runs by considering the corresponding probability dis- 
tributions in the next section. 

4.3.4. Structure of the interstellar medium 

The visual differences discussed above are quantified in 
Figure 15 where we show the cumulative mass fraction above 
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Fig. 14. — Face-on and edge-on maps of the galactic disk at t = 200 Myr showing gas surface density (top), mass weighted average gas temperature (middle) 
and stellar surface density (bottom). The face-on plots are calculated within z ± 1.5kpc of the disk to avoid excess halo material. Each panel is 24 kpc across. 
The temperature is calculated as f pT / f p, where the integral is performed along each pixel sightline within z± 150pc. The map of stellar distribution only 
includes star particles formed after the start of simulation, and does not include the star particles present in the initial conditions. 
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Fig. 15. — The cumulative mass fraction of the ISM above a given density n (left panel) and temperature T (right panel). The prominent tail at densities 
in excess of n > 100 cm -3 in the simulation without feedback is due to a population of dense, long-lived gas clumps. In the presence of feedback, such gas 
clumps are effectively dispersed which significantly reduces the fraction of gas at such densities. The temperature structure shows significant differences between 
different runs. In the run with delayed cooling significant « 10% of the disk gas is heated to temperature in excess of T ~ 10 5 K, while this fraction is only 
~ 0.1% in other runs. The "all" run has more mass around T ~ 10 3 K, which is associated with embedded star particles heating the ISM to warm temperatures. 
In the case of strong feedback, dense gas is not heated, but dispersed, and diffuse gas is heated to very high temperatures. This mechanism has little effect on the 
cumulative mass function, as the hot phase is negligible by mass. 




a given density and temperature at t = 200 Myr. All simula- 
tions are analyzed in the regions shown in Figure 14 within 
a distance of ±0.5 kpc of the disk plane. In the case of no 
feedback, the existence of dense gas clumps is manifested in 
the tail of the density distribution at ft > 100 cm -3 . The den- 
sity and temperature distributions in runs with feedback are 
qualitatively similar; the high-density tail at n w 100 cm -3 is 
suppressed as gas in star forming regions is efficiently dis- 
persed. Simulation with a second feedback energy variable 
has the least amount of dense gas, as could be deduced from 
its SFR in the bottom-right panel of Figure [TT] We note that 
the details of the high-density tail, as well astne the dispersal 
process, likely depend on the choice of star formation den- 
sity threshold and numerical resolution. The distributions pre- 
sented here are useful in interpreting trends of the KS relation 
normalization discussed above. For example, it is clear that 
runs with efficient feedback have SFR comparable to the run 



with no feedback and ten times lower because they simple 
have less dense gas. 

The temperature structure in the right panel also reveals sig- 
nificant differences between feedback schemes. In runs with 
delayed cooling, ~ 10% of the disk's gas mass is at T > 10 5 K, 
which is two orders of magnitudes greater than in the other 
runs. This can be seen in the temperature map in Figure 14 
where the central region features a hole of hot, ionized, but 
dense, gas formed out of percolating star forming regions of 
feedback ejecta. However, all runs have a comparable fraction 
of gas in the hot coronal phase (T > 10 6 K). For comparison, 
in the Milky Way dis k < 1 % of the gas mass is thought to be 
in the hot phase (e.g. Ferriere 2001 ). The prominent bump in 
the "all" run around T ~ 10 3 K is associated with embedded 
star particles heating the ISM to warm temperatures. In the 
strong feedback models "all_tau30" and "all_f05_tl0", feed- 
back disperses dense gas disperses more efficiently, and heat- 
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ing occurs in the diffuse rather then dense phase, which is 
why there is no significant mass contribution in the warm or 
hot phase from these runs in this figure. 

Figure shows the density, (dV/dlogft)/V to t, and temper- 
ature PDFs7 (dV/dlogr)/V tot , defined as a fraction of disk 
volume in a given density or temperature range. A log-normal 
PDF is not a good description to the density PDF in our sim- 
ulations contrary to results of Wada & Norman ( 2007}, al- 
though it may be possible to describe the PDFs as super- 
positions of several l og-normal distributions corre sponding to 
different gas phases (Robertson & Kravtsov 2008 ). The figure 
shows that the run without feedback has the most dense gas, 
but the smallest amount of tenuous gas at n < 0.1 cm -3 . Inter- 
estingly, the run with delayed cooling has less tenuous gas of 
density n ~ 10~ 2 - 10~ 3 cm -3 than our fiducial run. This indi- 
cates that feedback models with early feedback injection can 
efficiently create both a diffuse ionized warm phase and a ten- 
uous coronal phase without resorting to artificially delaying 
gas cooling. 

The multiphase structure of the ISM is apparent in the tem- 
perature PD F, where all simulations sho w signatures of a three 
phase ISM ( |McKee & Qstriker|[l977b| ), connected by gas at 
intermediate temperatures. Without feedback, the gas cools 
down to a very thin disk (only a few cells in vertical height) 
with a substantially lower contribution to the volume in the 
cold phase (T < 10 4 K) compare to runs with feedback, which 
all feature thicker cold gas disks due to feedback driven tur- 
bulence. In addition, more cold gas is lost in star formation 
events when feedback is absent. This discrepancy is espe- 
cially apparent when comparing to the most efficient feedback 
run, "all_f05_tl0". The hot (T ~ 10 6 K) tenuous gas phase is 
present in all runs, although vigorous heating in "all_dc40" 
and "all_f05_tl0" creates pockets of gas at ~ 10 7 K, which 
vent out of the disk to the surrounding corona. As can be seen 
in Figure [14| the circum-galactic medium is more structured 
in "all" ancPall_tau30", which is evident from the wider dis- 
tribution of gas at T ~ 10 4 - 10 8 K. 

4.3.5. Velocity dispersion profiles 

We quantify the level of turbulent gas motions in the disks 
via the mass weighted, vertical line-of-sight velocity disper- 
sion profile cr z (r), shown in Figure [T7| for the gas cold com- 
ponent (T < 10 4 K). Such profiles can be observed in real 
galaxies and comparisons of model results and observations 
can help to constrain parameters of feedback models. Indeed, 
we could expect that models with the most efficient feedback 
generate stronger gas motions, which should be manifested 
in larger velocity dispersions. The figure shows that signif- 
icant velocity dispersion declining with increasing radius is 
produced in all runs. Such declining dispersion profiles are 
indeed observed in spiral galaxies for the neutral HI gas (e.g. 
Meur er et aLl[T996l [Petric & Ruj^|2007| |Tamburro efaL] 
2009 ). The fact that significant velocity dispersion is observed 
in the run with no feedback, indicates that most of the mo- 
tions are due to disk instabilities and not due to feedback per 
se. In fact, velocity dispersion in the inner regions is even 
somewhat smaller in our fiducial run with the "All" feedback 
model. This difference is probably due to formation of mas- 
sive gas clumps which can more efficiently stir the gas as they 
move around and merge with each other in the weaker feed- 
back runs. Nevertheless, the largest velocity dispersions, in 
the inner 10 kpc of the disk, are observed in runs with de- 
layed cooling and large tir, i.e. models with the most efficient 
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Fig. 17. — The vertical line-of-sight velocity dispersion a z (r). The lines are 
described in Figure [15] All simulations show a radially declining dispersion 
profile, settling on cr z ~ 5kms _1 in the outer parts of the disk (r > 10 kpc) 
where Ssfr = 0. The central increase in turbulent dispersions occurs even in 
the case of no feedback (magenta line), illustrating the propensity of gravita- 
tional instabilities in generating random motions. As feedback is boosted, the 
velocity dispersions increase significantly towards the central, star forming, 
part of the galaxy. In the case of delayed cooling (red line), the dispersions 
are almost twice as high as in the standard "all" run (black line), and a similar 
trend, although slightly weaker, is found for tir = 30 (black dashed line), or 
a separate feedback energy variable is adopted (blue line). 

feedback. 

Using the THINGS galaxy sample, [Tamburro et al.| ( [20Q9 



analyzed the radial HI velocity dispersion, am, and star for- 
mation rate surface density profiles and found positive cor- 
relation between the kinetic energy of HI and the SFR. The 
increase in am at smaller radii indeed correlates with an in- 
crease in star formation activity, both in observations and 
simulations, but so does the level of shear and strength of 
disk self-gravity. Gravitational instabilities can generate a 
significant base line level o f turbulence even wi thout any 
contribution from feedback ( [Agertz et al.| [2009a), as illus 



trated in Figure [17 
plateau of am r 
(S SF r) < 10" 3 



Observations indicate a characteristic 
lOkms -1 in galaxies with a globally a veraged 
lO^Afoyr^kpc" 2 <Dib et al.||2006|), above 



which stellar feedback becomes the more dominant driver of 
the observed HI velocit y dispersions (as shown numerically 
by | Agertz et al.||2009a| ). The propensity for different feed- 
back models to generate turbulent velocity dispersions in ISM 
gas may therefore manifest more strongly in starbursting sys- 
tems. We leave an investigation of the velocity dispersion de- 
pendence on feedback parameters and star formation surface 
density for a future study. 

5. DISCUSSION AND CONCLUSIONS 

In this paper we have presented a new model for stellar 
feedback that explicitly considers the injection of both mo- 
mentum and energy in a time resolved fashion. In particular, 
we have calculated the time dependent momentum and energy 
budget from radiation pressure, stellar winds, supernovae type 
II and la, as well as the associated mass and metal loss for all 
relevant processes. We present a novel prescription for model- 
ing the early (pre-SNII) injection of momentum due to stellar 
winds and radiation pressure from massive young stars. These 
stellar feedback processes were implemented and tested in the 
AMR code RAMSES. We have also examined and compared 
the effects of feedback in this new implementation and other 
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popular recipes on properties of simulated galactic disks. 

Using idealized simulations of star forming patches of gas 
and star forming spiral galaxies, we study how each stellar 
feedback source affects the overall rate of galactic star forma- 
tion, as well as density, temperature, and velocity structure of 
the ISM. We find that early pre-SN injection of momentum is 
an important ingredient, which qualitatively changes the ef- 
fectiveness of stellar feedback. In a given stellar population, 
supernovae explode only after ~ 4Myr, while essentially all 
momentum and energy associated with radiation pressure and 
stellar winds are deposited in the first 3 - 4Myr. We show that 
such momentum injection disperses dense gas in star forming 
regions, which drastically increases the impact of subsequent 
SNII energy injection, even when no delay of cooling is as- 
sumed. Our simulations of massive (M ~ 1O 6 M ) star form- 
ing clouds indicate that momentum based feedback alone can 
limit the global cloud star formation efficiency to e c \ ~ 10%. 
In absence of the pre-SN momentum feedback, we re cover 



the classical over-cooling pro blem for stellar feedback ( Katz 
1992[|Navarro & W hite 1993 ), as gas cooling times are short 
in the dense star forming ISM (t coo \ ~ 10 3 years), and star for- 
mation is less affected by feedback. 

In a simulated Milky Way-like galaxy, we find that star for- 
mation rates, and the normalization of the Kennicutt- Schmidt 
relation, are significantly affected by inclusion of stellar feed- 
back. Interestingly, we find that the normalization of the 
Kennicutt-Schmidt relation is less sensitive to the assumed 
star formation efficiency per free-fall time (eff) in schemes 
with efficient feedback due to self-regulating effect of feed- 
back on density and temperature PDFs within interstellar 
medium of simulated galaxies. An order of magnitude change 
in eff only results in only a factor of two increase in the KS re- 
lation normalization. 

Our results illustrate the importance of not only accounting 
for the entire momentum and energy budget of stellar feed- 
back, but also to inject momentum and energy at the appro- 
priate stages of stellar evolution. A similar conclusion was 



recently reac hed by Hopkin s et al.| ( |20lTb| ) (see also Hopkins 
et al.||2012a| based on high-resolution SPH simulations. In 



this paper we show how this effect can be incorporated at the 
resolution typical for state-of-the-art galaxy formation simu- 
lations. 

Although the qualitative trends illustrated by our results are 
clear, it is not obvious whether the effects of feedback, es- 
pecially the survival and impact of shocked winds and SNe 
ejecta, are modelled correctly. This is because any subgrid 
feedback scheme by necessity is implemented at scales close 
to the resolution of the simulations, where numerical effects 
play a role. In our experiments we find that even if only 
~ 10% of thermal feed back energy is retained for 1 - lOMyr 
(as su ggested by e.g. Thornton et al.|[T998} |Cho & Kang 



2008), stored and followed using a separate energy variable, 
this energy has a significant effect on star formation rates, the 
ISM density structure and turbulent velocity dispersions. 

Comparing different feedback prescriptions, we find that 
the recipe presented in this paper results in effects on galactic 



star formation rate and interstellar medium structure similar 
to the results of feedback schemes with a delay of feedback 
energy dissipation if the infrared optical depth in star form- 
ing regions is sufficiently hi gh (ttr > 10). This con clusion is 
c onsistent with the result s of |Hopkins et al.| ( [201 lbj ). 

Hopkins et al. ( 2011b| reported average values of (ttr) ~ 
10-30 in their isolated "Milky Way" SPH simulation. This 
value refers to the typical ti R adopted as an SPH particle 
is kicked by their feedback scheme. In the empirically- 
motivated subgrid model for radiation pressure momentum 
presented in this paper, such high values of tir are achieved 
only around massive star clusters, M c \ > 1O 6 M , which are 
rare in our simulations of galactic disks. The actual values 
of IR optical depth around young clusters are quite uncertain. 
If large values, tir > 30, indeed are appropriate, this can be 
incorporated as a normalizati on con stant in our relation in Ap- 
pendix [A] (i.e. i]2 in Equation All| ). We note the that in some 
situations, even momentum from single scattering of photons 
(i.e. without IR trapping, tir ^ 0) can have a significant effect 



( [Wise et al |20T2| [Chattopadhya y'eTal.|2012| ). The regimes in 
which such feedback is efficient remain to be explored and 
clarified. 

The above conclusions are based on the simulations con- 
ducted at spatial resolutions typical of what is affordable by 
current cosmological simulations of galaxy formation, i.e. 
~ 10- lOOpc. We have not demonstrated numerical conver- 
gence in this work, and we do not necessarily expected this to 
be easily achieved; as resolution improves, the density PDF 
changes as self-gravitating gas can collapse to higher densities 
and the gas dissipates energy at a higher rate. This leads to a 
shorter star formation time scale (as tsF ~ P~ ' 5 , Equation 23 ) 
and hence an increased rate of star formation. Numerical con- 
vergence can in principle be achieved by imposing a pressure 
floor via a polytropic equation of state, P ex p 7 , where 7 = 2, 
similar to what is nec essary to avoid artificial fragmentation 
( Truelove et al.|[T997] ). In this case we impose, by hand, a 
floor to the allowed minimum Jeans mass, for which conver- 
gence in principle is achievable. However, in this case sim- 
ulations may converge to an incorrect (and arbitrary) result, 
if the polytropic equation of state does not capture the actual 
thermodynamic properties of ISM realistically. 

It is clear that any implementation of the star formation- 
feedback loop requires thorough testing against observations, 
such as the Kennicutt-Schmidt relation, velocity dispersion 
profiles of gas etc. We plan to carry out such tests using 
the implementations of the feedback models described in this 
paper, as well as different implementations of star formation 
recipes, in self-consistent cosmological galaxy formation sim- 
ulations in future work. 
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APPENDIX 

THE SUBGRID MODEL FOR RADIATION PRESSURE MOMENTUM 

One of the main difficulties in modelling momentum transferred to gas by radiation pressure is in proper accounting for 
contribu tion of momentum due to multiple scatterings of infrared photons by dust grains. In the implementation of Hopkins et al. 
( |2011a| ), an iterative clump finding algorithm was used to identify star forming clouds. All stars within the cloud radius transfer 
momentum to the gaseous components according to Equation [5] where the infrared optical depth r IR = ^ IR S gas . Here E gas is the 
gas surface density and Hopkins et al. adopt a constant opacity ft IR w 5 cm 2 g _1 , which is appropriate for dust temperatures of 
r d >100K. 

The surface density of gas in star forming clump was calculated directly from simulations as E c i = M c i/(7r^ 1 ), where R c \ is the 
radius given by the clump finding routine. Hopkins et al. report average values (at gas particle launch) of (tjr) ~ 10-30 in the 
Milky Way environment, and ~ 30- 100 in high redshift disk analogues. Values of this magnitude are a direct outcome of clouds 
collapsing to the point where densities are high enough for tir to halt the process, leading to cloud collapse. However, this process 
is highly sensitive to simulation resolution and other numerical effects. Moreover, surface density of real star forming clumps 
depend on internal processes within these regions, such as supersonic turbulence and feedback, which will not be resolved even 
with the ~ parsec resolution. This uncertainty hence propagates into the calculation of the momentum transfer via reprocessing 
of IR radiation by dust. 

Below we discuss an alternative way of estimating tir based on observed properties of young star clusters and molecular 
clumps, which can readily be implemented in simulations adopting spatial a resolution of Ax ~ 10- lOOpc. 



Surface densities of observed star forming clumps 

Observed molecular clouds have sizes of ~ 5 - 100 pc, average den sities of ~ 100 cm -3 , and surface densities of Egmc ~ 
50-200 M pc~ 2 (e.g., |Bolatto et al.|2008 |Fukui & Kawamura|2010| . They have complex internal structure with gas density 
I orders of magnitude. This structure is thought to arise due to supersonic turbule 



ranging over several orders of magnitude . This st ructure is thought to arise due to supersonic turbulent flows and gravitational 
contraction during cloud formation (e.g., |Padoan et al.| 19971 |Li et al.|2004||Kritsuk et al.|2007] ). 

Star clusters form in clumps which have densiti es of > 10 3 - 1 n4 ^™ 3 i^i^nd ctorc 

higher density within the clumps (> 1 5 cm -3 
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1994 



2000| [l 



Lada & Lada 



,e.g. 



10 4 cm 3 , while individual stars form in high-density cor es of even 



Lada & Lada 



2003 ). Clumps in the Milky Way molecular clouds (e.g., 



2003 ) have masses of up to few x 1O^M and their surface densities range from < 



Williams 
15OM pc" z 



to 

15000M o pc -2 (or ~ 0.03 - 3gcm~ 2 ). Based on the structural parameters of more massive young star clusters, their parent gas 
clumps had even higher surface densities. Young (t < lOMyr) sta r clusters Milky Wa y clusters with masses of order ~ 1O 4 M 
are have half light radii of r^ ~ 0.5-2 pc (see compilation by Portegies Zw art et aT]|2010|). Observations of young clusters in 
external galaxies find r^ ~ 2pc independent of luminosity or cluster mass (M51, Scheepmaker et al. 2007 ) and young super star 
clusters in star burst galaxies e.g. M82 show similar properties ( [McCrady & Gr aham 2007 a|). 

The bottom panel of Figure [18] shows the relation between half -mass radius and mass for molecular clumps in the Milky Way 
(from the compilation shown in Figure 1 of [Fall et al.|2010b|» and for yo ung (age < 2 x 10 7 yrs) star clusters in the Milky Way 
and other nearby galaxies from the compilation of |Portegies Zwart et al.| ( [201 0| > as a function of their mass. Different lines show 
relations derived for star clusters and clumps in several recent studies, as described in the figure caption. The corresponding 
surface densities of the clumps and star clusters are plotted in the upper panel, and shows that although scatter is substantial, 
clumps and clusters in the MW follow a similar relation at masses < 10 5 M : R ex M a with a ~ 0.3-0.5. This implies that it is 
reasonable to assume that radii and masses of young clusters are a good reflection of the corresponding properties of their parent 
molecular clumps. 

For clusters of mass > 10 5 M the relation flattens (a w 0), although the scatter is large. The broken dashed line shows an 
approximation to the observed behavior of clumps and clusters: 



Rd=( 



V 3000M^ 



0.4 



pc, for M cl < 3 x 1O 4 M , 



tf cl = 2.5 pc, for M cl > 3 x 1O 4 M . 



(Al) 
(A2) 



The corresponding data and lines for surface densities defined as E d = M c \/ (27r7?^) are shown in the upper panel of the figure. For 
M c i < 3 x 10 4 Mq the surface densities are generally E < 1 gem -2 , while for more massive clusters they can reach E > 10 gem -2 . 
For the latter values of surface densities, the optical depth tir > 50 if dust temperatures are warm (T d > 200 K), as illustrated in 
Figure [3] 
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log M (M ) 



Fig. 18. — Surface dens ity versus mass (top panel) and half-mass radius versus mass (bottom panel) for the m olecular clumps in the Milky Way (blue squares, 
from the compilation by Fall et al. 2010b in their Figure 1) and young star clusters (stars, from compilation by Portegies Zwart et al. 2010). The green stars show 
clusters in the LMC ( Mackey & Gilmore 2003) and SMC, while cyan stars show star clu sters within Milky Way. The other star symbols show clusters in other 
galaxies, including starbursts such as M82 (solid red stars, Krumholz & Matzner 2009b) and the Antennae galaxies. Note that surface densitie s are estimated 
within the half mass radius: S = M/(2ttR 2 ). The broken magenta dashed lines show power law approximation to the clusters given by Equation |A2| Solid line 
in the bottom pa nel is power law fit to the MW low-mass star clusters from Lada & Lada (2003 1, while dotted line is fit to the mass-radius of Mw clumps from 
|Dib et al.| ( [2010) . 

One complication to considering radii and masses of observed young clusters is that they can evolve with time from the time 
of their birth, when most of the radiative pressure feedback has occurred. Indeed, observations show some weak correlation of 
cluster sizes with age (e.g., see Portegies Zwart et aL]|2010| for a recent review; in particular their Figure 8). Figure [T9| shows 
the masses a nd half-mass radii of young star clusters m the Local Group and beyond, from the compilation of Portegies Zwart 
et al.| ( |2010| ), as a function of their age. While some correlation of radii with age is apparent, the mass-age panel shows that 



interpretation of this correlation is not straightforward. While youngest clusters do have somewhat smaller ages, they also have 
smallest masses. 

Interestingly, the figure shows that there are quite a few observed massive clusters (M c \ > 10 5 M ) with ages ~ 5-20 x 10 6 yrs. 
Such clusters are largely missing, however, at smaller ages. This is likely because most of the distant clusters of such mass are 
deeply embedded in gas and dust and therefore have been misse d in observations. Some of the massive young stellar clusters are 
indeed observed to be deeply embedded in dust (A v ~ 9 - 10 mag | Gilbert et al.|2000[|Gilbert & G raham 2007 ). The most massive 
cluster in the Antennae galaxies, for example, is very faint in the optical band, but is one of the brighest sources in the IR. This 
indicates that conditions for substantial radiation pressure do exist in the natal gas clumps. 

Another interesting feature of the figure is the fact that there are almost no massive clusters with ages > 20 Myr. Unless there 
is some selection effect, this probably indicates that massive clusters dissolve as a result of dispersal of their parent clump gas 
and subsequent stellar mass loss. 

In summary, it is not obvious that structural parameters of massive clusters at age ~ 10 7 yrs are significantly different from the 



parameters of these clusters at birth. It is thus premature to apply an age correction to the data shown in the Figure 1 8 



A sub grid model for radiation pressure momentum 

As the above discussion illustrates, a direct implementation of radiation pressure momentum injection is not feasible in galaxy 
formation simulations, as it requires a resolved density structure of star forming clouds at parsec scales. For more general 
application, it is useful to develop a subgrid model based on empirical knowledge of structure and physics of star forming 
molecular clumps, which could be valid for different spatial resolutions. The proposed model is local in it is nature, acting only 
in the local resolution elements surrounding young stars, and therefore does not account for continuous acceleration of gas far 



outsid e of galactic disks, as recently proposed by Murray et al. (201 1 ), and modelled in the radiative transfer simulations by Wise 
|eTS1 ( |20T2> , 

In this model, any star particle formed by star formation recipe in simulations is regarded as an ensemble of star clusters, with 
an associated ensemble of natal molecular clumps onto which radiation pressure acts at early times. Such model can then be 
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Fig . 19. — Mass (top panel) and half-mass radius (bottom panel) of young star clusters as a function of their age using cluster compilation of Portegies Zwart 
|et al.|{20ToV The blue stars show clusters with ages longer than three dynamical times, while cyan stars show clusters with ages shorter than three dynamical 
times. The latter are classified as stellar associations. 



used to calculate the total momentum input from stars in all clumps within star forming cell by integrating over the clump mass 
function. 

We start off by defining the rate of momentum deposition imparted by radiation pressure from young stars on a molecular 
clump as 

Lit) 

Pel = (Vi + mrm) , (A3) 

c 

where L(t) is the bolometric luminosity of stars in a star cluster of age t, and as before ttr = ^irS c i is the IR optical depth. 
The first term describes the direct radiation absorption/scattering and should in principle be oc [1 -exp(-ruv)L but given that 
UV optical depth is always very large in dense star forming regions, r]i w 1. The optical and UV photons heat dust particles in 
surrounding gas and IR photons radiated by dust can transfer additional momentum if gas is optically thick to the IR radiation. 
The second term thus sp ecifies momentum transferred via multiple scatterings of IR photons re-radiated by dust particles (see, 
e.g. |Gayley et al.||1995] ), where 772 is added to parametrize possible modifications to the adopted ttr. As for 771, our fiducial 
choice is r?2 = 1, although factors of a few maybe be motivated due to grid smearing/cancelations for large advection velocities, 
as discussed below. 

By integrating p c \ over the clump mass function, we obtain the total imparted momentum rate from all star clusters onto their 
natal clumps, 



Pxox 



^ Mcl^min 



p c \i/j(M c i)dM ch 



(A4) 



where the minimum and maximum clump masses (M c i m i n and M c i >max ) set the normalization of the cluster mass function. We 
approximate the observed mass function of molecular clumps by a power law 



^(M cl )=A cl Mj 



(A5) 



with /? « 1.7 ±0.2 dKramer et al.||1998] ) similar to the power law slope of the molecular clouds themselves (e.g., Fukui & 
Kawamura 2010]). The latter can be approximated by a Schechter li ke function with exponential cutoff. The mass function of 
young star clusters can also be approximated by the Schechter form ([Zhang & Fall|19 99 ; de Grijs et al.|2003[ |Bik et al.]2 003 1 
Cresci etaL]|2005| |McCrady & Gra ham 2007b] see |Portegies Zwart et al.| |2010| for a review). Different studies that sample 
different parts of cluster mass function and often approximate it with a simple power law, can get somewhat different values of 
the slope. Nevertheless, the mass function of star clusters is generally found to be quite similar in shape to that of the molecular 
clouds and clumps. The similarity of mass function slopes implies that star formation efficiency in clumps, e c \ = M* jC i/M c i 
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(where M* c \ is mass of a star cluster a given clump of mass M c \ forms), is approximately independent of clump mass (Fall et al. 
|2010b| ). 

We interpret a formed star particle of mass m* as the total mass of an ensemble of star clusters formed with the constant 



efficiency e c \ from molecular clumps with a mass function given by Equation |A5) Although the actual mass function may have 
Schechter form, for our purposes we can approximate it as a power law (i.e., small-mass end of the Schechter function) with 
effective maximum clump mass, M c i max . The choice of this maximum mass is related to the normalization of the cluster mass 
function and we will assume that the maximum mass of star cluster M* iC \ imax = e c i^ci,max is 

^*,cl,max = Mmax^*- (A6) 

Motivation for such relation comes from observation that masses of young (ages of < 10 7 yrs) star clusters in M33 and LMC 



are ~ 10-50 times smaller than the masses of their parent molecular clouds (Fu kui & Kawamura|2010[ Portegies Zwart et al. 



2010| ). Given that m* is related to the mass of molecular gas in the cell via the star formation recipe, Equation |A6| establishes a 



relation between the masses of largest star cluster and mass of the molecular gas in the cell. Depending on assumed star formation 



efficiency, the value /i max ~ 0.1 - 1 is reasonable. The minimum mass of the clumps can be taken to be Af R i,min ~ 100 M (Lada 
& Lada 2003). In this work we fix the slope of the mass function to the above suggested value of f3 = 1.7. For this slope, most or 



the mass of clumps is in most massive clumps. This will maximize radiation pressure feedback compared to mass functions that 
have slopes fj > 2. 

For a clump of mass M c \ and half-mass radius R^, we can define gas surface density E c i = (M/2)/(ttR^) 9 velocity dispersion 
a = escape velocity v Q = 2a, and crossing time t c = R^/a, all of which could be fully characterized in terms of 

mass if we assume that clump mass and radius are related via power law: 

R h = C R M%. (Al) 

The slope of this relation is constrained by observed slope of the average E c i- M c \ relation, 1 -2a, and corresponding relation 



for young star clusters as indicated in Figure 1 8 In this work we adopt the relation given by Equation A2 



We obtain the luminosity L(t) from STARBURST99, and define 

L(t) = L x (t)M^ ch (A8) 
where L\(t) is bolometric luminosity per M . L\ is approximately constant at w 3.8 x 10 36 ergss -1 Mq 1 for t < 3 x 10 6 yrs, and 



A4 



decreases roughly as ~ t~ 125 at late times, see Figure [I] Using the above relations, the total momentum rate p tot in Equation 
may now be evaluated. The relation has two terms; let us consider them in turn. 

The first term is independent of surface density and, when integrated over the clump mass function, will simply give 

Aot,i =m m *; (A9) 

c 

where t is the age of a given stellar particle of mass m* in cell under consideration. This contribution can simply be summed up 
for all young stellar particles in the cell with ages as old as it is deemed to be significant, typically for a few Myr. 

The second term depends on the surface density. Before the parent molecular clumps are dispersed by their child star clusters, 
i.e. for time less than clump lifetime t < t c \, radiation pressure operates on the surface density of the clumps, 

S d = (l-e d )^%, (A10) 

where (1 - e c \) factor takes into account the fact that fraction e c \ of clump mass was turned into stars, while factor of 0.5 takes 
into account that R^ in the assumed mass-radius relation is the half-mass radius. Thus, for t < t c \, using the equations above, and 
integrating over the clump mass function, we obtain 

. 7/2ttlR(l-gcl)(2-/3) / /imax V" 2 " 1 ~ (M clM J M cl ^ x ) 3 ' 2 ^ L X (t) 2(l _ a) 

Aot ' 2a 2kCI 3 -2a- (3 \ e d J 1 -(M cl , mm /M cl , max )^ c m * ' ^ Ai ^ 

For t > t c \, the clump is dispersed and the radiation pressure will simply act on the surface density of the cell, Hg as ,c, with a 
possible boosting by some clumping factor to account for a clumpy nature of parent molecular cloud. The latter can be introduced 
via 772. The total momentum rate in this case will therefore be: 

Aot,2b = ^2^IRS g as,c^* (A12) 
C 

Summarizing, the total momentum rate is 

• t _ f Aot,l +Aot,2a if ^ < t c \, (A13) 

t0t 1 Aot,l + Aot,2b = (rjl + ^2^IRS gas ,c)^* ^ if t > t ch 

The clump lifetime t c \ is a highl y uncertain factor, b ut can reasonably be assumed to be a fixed multiple of the clump crossing 
time: t c \ = N c t c with N c ~ 5 - 10 (Palla & Stahler 2000), where crossing time t c is given by 

3 /2 

Rh_ C R (a-l)/2 , . . „ 



26 



Agertz et al. 



where R = CrM% relation was used. For massive clusters, Cr & 2.5 pc and a w and 

^ = 9 - 3xl ° 4 (i^) yrs - (A15) 

Thus, for N c ~ 5 - 10, th e life time of a clump is 0.5 - 1 x 10 6 yrs. However, the first stage of cluster life before clump dispersal 
is highly uncertain (e.g., Porte gies Zwart et al.|2010| ) and we do not know neither the exact lifetime nor its scaling with cluster 
mass. 

In the current work, we adopt the above model for radiation pressure feedback assuming t c \ = 3Myr, e c \ = 0.2, /i max = 1, 
ftiR = 5 (Z /Z) cm 2 g -1 and 771 = 772 = 2 for our fiducial model, unless noted otherwise. The metallicity scaling on the opacity 
is a crude way of accounting for the varying dust-to-gas ratios. The values of 77 account for the fact that the actual measured 
momentum injected by a star particle is found to be reduced during advection of gas through computational mestf^] 

Caveats related to the stellar mass used for radiation pressure 



The particle mass m* entering the terms of Equation | A 1 3 1 does not necessarily need to be interpreted as the mass of each star 
particle. In fact, this is not preferred as the strength of radiation pressure depends on m* in a non-linear fashion. As the particle 
mass usually is a function or numerical resolution, and the effect of radiation pressure would weaken at higher resolution when 
the numerical scheme allows for the formation of lower mass star particles. When calculating the effect of radiation pressure we 
therefore adopt 

n 

rn* — » m* = ^m* ?n (0 for t <t ch (A16) 
i=i 

i.e. m* is the binned mass of all n star particles in a cell younger than some given age. As indicated in the above relation,we 
simply set this age to be the same as the clump life time t c \ defined above, for which we adopt the fiducial value of 3 Myr. This is 
consistent with estimates of the duration of embedded stage of young clusters (e.g., |Portegies Zwart et al.|2010| ). 

Another caveat in our implementation of radiation pressure is that we only consider the young stars in one computational cell 
when estimating r IR . The actual cell has nothing to do with the physics of the problem, and is just a convenient implementation 
choice. At very high resolution, Ax ~ few pc, star formation will be spread over several cells in massive GMCs and the collective 
effect of radiation pressure will be underesti mated. This issue can b e avoided by searching neighboring cell for young stars, or 
by using a star cluster finder, as suggested by Hop kins et al.| ( |20lTb| ). In the regime of Ax ~ 50- lOOpc, adopted for the galactic 
disks in § |4.3| where cell sizes matches observed sizes of massive GMCs, this is not an issue. 

IMPLEMENTATION OF NON-THERMAL PRESSURE 

As discussed in §[3] feedback momentum can be injected directly to the computational grid via "kicks." An alternative approach 
is to allow for the hydro scheme to generate momentum by appropriately pressurizing the local volume in which momentum is 
injected. In this approach, which we adopt in a sub-set of our simulations, we define the non-thermal pressure due to injected 
momentum as 

Pnt=P/A (Bl) 

where the area A is an arbitrary computational region, here chosen to be the surface area of computational cell (A = 6 Ax 2 ) 
containing a young star particle used to compute the momentum injection in the subgrid model described above. This pressure is 
added to the effective pressure, 

P&ff = P therm +^nt ? (B2) 

where P t herm is the thermal pressure. P e ff replaces P t herm in the sound speed definition when calculating the time step, and is 
otherwise only actively involved in the flux calculation (i.e., in the Godunov step). Here P e ff is consistently traced to the cell 
interfaces (here using the piecewise linear approximation) as a separate pressure variable, using its own TVD slopes. The left and 
right-hand states of P e ff then enter the Riemann solver, where it replaces the thermal pressure. Specifically, the momentum and 
energy equations that we aim to solve are 

^-(pv) + V-(pv® V+P&) = -pV(j> (B3) 
ot 

and 

d 

—(pE) + V- [pv(E+P Qff /p)] = -pv • V0, (B4) 

where E is the specific total energy, <\> the gravitational potential, and ® is the outer vector product. The effective pressure hence 
never enters in to the specific total energy (as P e ff/p(7~_l)) w hich is an important distinction to make in the MUSCL-scheme 
( |van Leer|1979| ) adopted in the RAMSES code (Teyssier 2002]), as it traditionally only stores one variable representing the total 



energy (the conservative variable), and pressure (the primitive variable) is derived from it after subtracting the kinetic energy. To 

15 This is due to two effects: smearing occurring at the grid level and be Qn ^ Qrder Qf „ 15 _ 25% M dimension for translationa i ve i oc i- 

momentum cancellations introduced by the nearest grid point approach. The . _ x 

latter occurs as a star particle discretely switches cells and reduce the momen- ties ( relatlve t0 tne g nd ) U P t0 v ~ lOOUkms . 
turn flux injected from a previous time-step. We have measured this effect to 
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universally replace P t herm by P e ff everywhere in the method would hence not be consistent with the equations we want to solve 
(as well as in the cooling routines). For our implementation of radiation pressure, the effective pressure is not advected, but is 
updated every fine time-step in the feedback routine, and stored as a se parat e variable. 

In the case of the second feedback energy variabl e Efp , introduced in § 3.2 we calculate a non-thermal pressure P nt = (7- l)p£fb, 
which enters the effective pressure as in Equation |B2| This quantity is then treated exactly as described above, although Er> is 
passively advected with the flow, i.e. it obeys 

^(pE fb ) + V.(pvE fb ) = 0. (B5) 



